hw04 done

.gitignore

@@ -1,2 +1,3 @@
 build
 GNUmakefile
+.vscode

.vscode/settings.json

@@ -0,0 +1,45 @@
+{
+    "files.associations": {
+        "cmath": "cpp",
+        "deque": "cpp",
+        "string": "cpp",
+        "vector": "cpp",
+        "array": "cpp",
+        "atomic": "cpp",
+        "*.tcc": "cpp",
+        "cctype": "cpp",
+        "chrono": "cpp",
+        "clocale": "cpp",
+        "cstdarg": "cpp",
+        "cstddef": "cpp",
+        "cstdint": "cpp",
+        "cstdio": "cpp",
+        "cstdlib": "cpp",
+        "ctime": "cpp",
+        "cwchar": "cpp",
+        "cwctype": "cpp",
+        "unordered_map": "cpp",
+        "exception": "cpp",
+        "algorithm": "cpp",
+        "memory": "cpp",
+        "memory_resource": "cpp",
+        "optional": "cpp",
+        "ratio": "cpp",
+        "string_view": "cpp",
+        "system_error": "cpp",
+        "tuple": "cpp",
+        "type_traits": "cpp",
+        "utility": "cpp",
+        "fstream": "cpp",
+        "initializer_list": "cpp",
+        "iosfwd": "cpp",
+        "istream": "cpp",
+        "limits": "cpp",
+        "new": "cpp",
+        "ostream": "cpp",
+        "sstream": "cpp",
+        "stdexcept": "cpp",
+        "streambuf": "cpp",
+        "typeinfo": "cpp"
+    }
+}

CMakeLists.txt

@@ -2,8 +2,7 @@ cmake_minimum_required(VERSION 3.12)
 project(hellocmake LANGUAGES CXX)
 
 set(CMAKE_CXX_STANDARD 17)
-if (NOT CMAKE_BUILD_TYPE)
-    set(CMAKE_BUILD_TYPE Release)
-endif()
+set(CMAKE_BUILD_TYPE Release)
 
 add_executable(main main.cpp)
+target_compile_options(main PUBLIC -ffast-math -march=native)

main.cpp

@@ -4,66 +4,83 @@
 #include <chrono>
 #include <cmath>
 
+int randMax = 1 / RAND_MAX;
+
+const size_t N = 48;
+
 float frand() {
-    return (float)rand() / RAND_MAX * 2 - 1;
+    return (float)rand() *randMax * 2 - 1;
 }
 
 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 stars;
+
 
 void init() {
+    float tmp = frand();
+    #pragma GCC unroll 16
     for (int i = 0; i < 48; i++) {
-        stars.push_back({
-            frand(), frand(), frand(),
-            frand(), frand(), frand(),
-            frand() + 1,
-        });
+        stars.px[i] = tmp;
+        stars.py[i] = tmp;
+        stars.pz[i] = tmp;
+        stars.vx[i] = tmp;
+        stars.vy[i] = tmp;
+        stars.vz[i] = tmp;
+        stars.mass[i] = tmp+1;
     }
 }
 
 float G = 0.001;
 float eps = 0.001;
 float dt = 0.01;
+float Gdt = G * dt;
+float eps2 = eps * eps;
 
 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 < 48; i++) {
+        #pragma GCC unroll 16
+        for(size_t j = 0; j < 48; 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 + eps2;
+            d2 *= std::sqrt(d2);
+            float dao = 1 / d2;
+
+            stars.vx[i] += dx * stars.mass[j] * (Gdt * dao);
+            stars.vy[i] += dy * stars.mass[j] * (Gdt * dao);
+            stars.vz[i] += dz * stars.mass[j] * (Gdt * dao);
         }
-    }
-    for (auto &star: stars) {
-        star.px += star.vx * dt;
-        star.py += star.vy * dt;
-        star.pz += star.vz * dt;
+        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;
-            float d2 = dx * dx + dy * dy + dz * dz + eps * eps;
-            energy -= other.mass * star.mass * G / sqrt(d2) / 2;
+    Star other = stars;
+    Star star = stars;  
+    for (std::size_t i = 0;i <48;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 *0.5f;
+        #pragma GCC unroll 16
+        for (std::size_t j = 0;j < 48;j++) {
+            float dx = other.px[j] - star.px[i];
+            float dy = other.py[j] - star.py[i];
+            float dz = other.pz[j] - star.pz[i];
+            float d2 = dx * dx + dy * dy + dz * dz + eps2;
+            float tmp = 1 / sqrt(d2);
+            energy -= other.mass[j] * star.mass[i] * G * tmp *0.5f;
         }
     }
     return energy;
+
 }
 
 template <class Func>