finished hw04

CMakeLists.txt

@@ -7,3 +7,9 @@ if (NOT CMAKE_BUILD_TYPE)
 endif()
 
 add_executable(main main.cpp)
+
+if (MSVC) 
+    target_compile_options(main PUBLIC /fp:fast /arch:AVX2)
+else()
+    target_compile_options(main PUBLIC -ffast-math -march=native)
+endif()

main.cpp

@@ -4,63 +4,73 @@
 #include <chrono>
 #include <cmath>
 
+constexpr size_t star_size = 48;
+constexpr float G = 0.001;
+constexpr float eps = 0.001;
+constexpr float dt = 0.01;
+constexpr float eps2 = eps * eps;
+constexpr float Gdt = G * dt;
+
 float frand() {
     return (float)rand() / RAND_MAX * 2 - 1;
 }
 
+template<size_t N>
 struct Star {
-    float px, py, pz;
-    float vx, vy, vz;
-    float mass;
+    float px[N], py[N], pz[N];
+    float vx[N], vy[N], vz[N];
+    float mass[N];
 };
 
-std::vector<Star> stars;
+Star<star_size> star;
 
 void init() {
-    for (int i = 0; i < 48; i++) {
-        stars.push_back({
-            frand(), frand(), frand(),
-            frand(), frand(), frand(),
-            frand() + 1,
-        });
+    for (int i = 0; i < star_size; i++) {
+        star.px[i] = frand();
+        star.py[i] = frand();
+        star.pz[i] = frand();
+        star.vx[i] = frand();
+        star.vy[i] = frand();
+        star.vz[i] = frand();
+        star.mass[i] = frand() + 1;
     }
 }
 
-float G = 0.001;
-float eps = 0.001;
-float dt = 0.01;
-
 void step() {
-    for (auto &star: stars) {
-        for (auto &other: stars) {
-            float dx = other.px - star.px;
-            float dy = other.py - star.py;
-            float dz = other.pz - star.pz;
-            float d2 = dx * dx + dy * dy + dz * dz + eps * eps;
-            d2 *= sqrt(d2);
-            star.vx += dx * other.mass * G * dt / d2;
-            star.vy += dy * other.mass * G * dt / d2;
-            star.vz += dz * other.mass * G * dt / d2;
+    for (size_t i = 0; i < star_size; ++i) {
+        float sumx = 0.0f, sumy = 0.0f, sumz = 0.0f;
+        for (size_t j = 0; j < star_size; ++j) {
+            float dx = star.px[j] - star.px[i];
+            float dy = star.py[j] - star.py[i];
+            float dz = star.pz[j] - star.pz[i];
+            float d2 = dx * dx + dy * dy + dz * dz + eps2;
+            d2 *= std::sqrt(d2);
+            sumx += dx * star.mass[j] * Gdt / d2;
+            sumy += dy * star.mass[j] * Gdt / d2;
+            sumz += dz * star.mass[j] * Gdt / d2;
         }
+        star.vx[i] += sumx;
+        star.vy[i] += sumy;
+        star.vz[i] += sumz;
     }
-    for (auto &star: stars) {
-        star.px += star.vx * dt;
-        star.py += star.vy * dt;
-        star.pz += star.vz * dt;
+    for (size_t i = 0; i < star_size; ++i) {
+        star.px[i] += star.vx[i] * dt;
+        star.py[i] += star.vy[i] * dt;
+        star.pz[i] += star.vz[i] * dt;
     }
 }
 
 float calc() {
     float energy = 0;
-    for (auto &star: stars) {
-        float v2 = star.vx * star.vx + star.vy * star.vy + star.vz * star.vz;
-        energy += star.mass * v2 / 2;
-        for (auto &other: stars) {
-            float dx = other.px - star.px;
-            float dy = other.py - star.py;
-            float dz = other.pz - star.pz;
-            float d2 = dx * dx + dy * dy + dz * dz + eps * eps;
-            energy -= other.mass * star.mass * G / sqrt(d2) / 2;
+    for (size_t i = 0; i < star_size; ++i) {
+        float v2 = star.vx[i] * star.vx[i] + star.vy[i] * star.vy[i] + star.vz[i] * star.vz[i];
+        energy += star.mass[i] * v2 / 2;
+        for (size_t j = 0; j < star_size; ++j) {
+            float dx = star.px[j] - star.px[i];
+            float dy = star.py[j] - star.py[i];
+            float dz = star.pz[j] - star.pz[i];
+            float d2 = dx * dx + dy * dy + dz * dz + eps2;
+            energy -= star.mass[j] * star.mass[i] * G / std::sqrt(d2) / 2;
         }
     }
     return energy;