From 70bf32f9d480e771b1227d6cfc045b92841c9ec1 Mon Sep 17 00:00:00 2001
From: QifanWang <one302755067@gmail.com>
Date: Sat, 26 Nov 2022 20:20:45 +0800
Subject: [PATCH] finish hw04

---
 CMakeLists.txt |  1 +
 main.cpp       | 96 +++++++++++++++++++++++++++++++-------------------
 2 files changed, 61 insertions(+), 36 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 29b152c..6c85c06 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,3 +7,4 @@ if (NOT CMAKE_BUILD_TYPE)
 endif()
 
 add_executable(main main.cpp)
+target_compile_options(main PUBLIC -ffast-math -march=native)
diff --git a/main.cpp b/main.cpp
index cf6369b..7482e41 100644
--- a/main.cpp
+++ b/main.cpp
@@ -8,59 +8,83 @@ float frand() {
     return (float)rand() / RAND_MAX * 2 - 1;
 }
 
-struct Star {
-    float px, py, pz;
-    float vx, vy, vz;
-    float mass;
-};
+const int MAXLEN = 48;
 
-std::vector<Star> stars;
+struct Star {
+    float px[MAXLEN], py[MAXLEN], pz[MAXLEN];
+    float vx[MAXLEN], vy[MAXLEN], vz[MAXLEN];
+    float mass[MAXLEN];
+} stars;
 
 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 < MAXLEN; i++) {
+        stars.px[i] = frand();
+        stars.py[i] = frand();
+        stars.pz[i] = frand();
+        stars.vx[i] = frand();
+        stars.vy[i] = frand();
+        stars.vz[i] = frand();
+        stars.mass[i] = frand() + 1;
     }
 }
 
-float G = 0.001;
-float eps = 0.001;
-float dt = 0.01;
+const float G = 0.001;
+const float eps = 0.001;
+const float dt = 0.01;
+
+const float Gdt = G * dt;
 
 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;
+    for (size_t i = 0; i < MAXLEN; ++i)
+    {
+        for (size_t j = 0; j < MAXLEN; ++j)
+        {
+            float dx = stars.px[j] - stars.px[i];
+            float dy = stars.py[j] - stars.py[i];
+            float dz = stars.pz[j] - stars.pz[i];
             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;
+            d2 = 1 / std::sqrt(d2);
+            d2 = d2 * d2 * d2;
+            stars.vx[i] += dx * stars.mass[j] * Gdt * d2;
+            stars.vy[i] += dy * stars.mass[j] * Gdt * d2;
+            stars.vz[i] += dz * stars.mass[j] * Gdt * d2;
         }
     }
-    for (auto &star: stars) {
-        star.px += star.vx * dt;
-        star.py += star.vy * dt;
-        star.pz += star.vz * dt;
+    
+    // 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 = 1 / std::sqrt(d2);
+    //         d2 = d2 * d2 * d2;
+            
+    //     }
+    // }
+
+    for(size_t i = 0; i < MAXLEN; ++i) {
+        stars.px[i] += stars.vx[i] * dt;
+        stars.py[i] += stars.vy[i] * dt;
+        stars.pz[i] += stars.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;
+    for (size_t i = 0; i < MAXLEN; ++i)
+    {
+        float starvx = stars.vx[i], starvy = stars.vy[i], starvz = stars.vz[i];
+        float starpx = stars.px[i], starpy = stars.py[i], starpz = stars.pz[i];
+        float starmass = stars.mass[i];
+        float v2 = starvx * starvx + starvy * starvy + starvz * starvz;
+        energy += starmass * v2 / 2;
+        for (size_t j = 0; j < MAXLEN; ++j) {
+            float dx = stars.px[j] - starpx;
+            float dy = stars.py[j] - starpy;
+            float dz = stars.pz[j] - starpz;
             float d2 = dx * dx + dy * dy + dz * dz + eps * eps;
-            energy -= other.mass * star.mass * G / sqrt(d2) / 2;
+            energy -= stars.mass[j] * starmass * G / sqrt(d2) / 2;
         }
     }
     return energy;