diff --git a/Loop1/src/main.cpp b/Loop1/src/main.cpp index 318f940..af0498d 100644 --- a/Loop1/src/main.cpp +++ b/Loop1/src/main.cpp @@ -4,98 +4,158 @@ #include #include #include +#include -void calc(double* arr, uint32_t ySize, uint32_t xSize, int rank, int size) -{ - if (rank == 0 && size > 0) - { - for (uint32_t y = 0; y < ySize; y++) - { - for (uint32_t x = 0; x < xSize; x++) - { - arr[y*xSize + x] = sin(0.00001*arr[y*xSize + x]); - } - } - } -} +#include +#include -int main(int argc, char** argv) -{ - int rank = 0, size = 0, buf = 0; - uint32_t ySize = 0, xSize = 0; - double* arr = 0; - - MPI_Init(&argc, &argv); - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &size); - - if (rank == 0) - { - // Check arguments - if (argc != 3) - { - std::cout << "[Error] Usage \n"; - buf = 1; - MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); - return 1; - } +#define ROOT 0 - // Prepare input file - std::ifstream input(argv[1]); - if (!input.is_open()) - { - std::cout << "[Error] Can't open " << argv[1] << " for write\n"; - buf = 1; - MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); - return 1; +#define S_TAG 1 + +#define coord(x, y, xSize) ((y)*(xSize) + x) + +// arr[y*xSize + x] = sin(0.00001*arr[y*xSize + x]); +void calc(double* arr, uint32_t ySize, uint32_t xSize, int rank, int size) { + int32_t step = 0; + int32_t* range = NULL; + int32_t* displace = NULL; + + MPI_Bcast(&xSize, 1, MPI_UNSIGNED, ROOT, MPI_COMM_WORLD); + MPI_Bcast(&ySize, 1, MPI_UNSIGNED, ROOT, MPI_COMM_WORLD); + + if(rank == ROOT) { + uint32_t y_step = ceil(ySize / (1.0 * size)); + uint32_t prev_range = 0; + + range = (int32_t*)calloc(size, sizeof(int32_t)); + displace = (int32_t*)calloc(size, sizeof(int32_t)); + + for(int send_rank = 0; send_rank < size; ++send_rank) { + // We will send start_y, step, xSize and part of arr + uint32_t real_step = std::min(y_step, ySize - prev_range); + range[send_rank] = real_step * xSize; + displace[send_rank] = prev_range * xSize; + prev_range += y_step; + } } + MPI_Scatter(range, 1, MPI_INT, + &step, 1, MPI_INT, + ROOT, MPI_COMM_WORLD); - // Read arguments from input - input >> ySize >> xSize; - MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); - arr = new double[ySize * xSize]; + double* my_copy = (double*) calloc(step, sizeof(double)); + MPI_Scatterv(arr, range, displace, MPI_DOUBLE, + my_copy, step, MPI_DOUBLE, + ROOT, MPI_COMM_WORLD); - for (uint32_t y = 0; y < ySize; y++) - { - for (uint32_t x = 0; x < xSize; x++) - { - input >> arr[y*xSize + x]; - } + + for(int32_t x = 0; x < step; ++x) { + my_copy[x] = sin(0.00001*my_copy[x]); } - input.close(); - } else { - MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); - if (buf != 0) - { - return 1; + /* + if(size > 1 && rank == 1) { + for(int i = 0; i < step; ++i) { + printf("%lf \n", my_copy[i]); + } + }*/ + /* + if(rank == ROOT) { + for(int i = 0; i < size; ++i) { + printf("rang[i] = %d, dit[i] = %d\n", range[i], displace[i]); } - } + }*/ + + //printf("[%d] xSize = %u, step = %u, my_copy = %p\n", rank, xSize, step, my_copy); + + MPI_Gatherv(my_copy, step, MPI_DOUBLE, + arr, range, displace, MPI_DOUBLE, + ROOT, MPI_COMM_WORLD); + + if(rank == ROOT) { + + free(range); + free(displace); + } + free(my_copy); +} + +int main(int argc, char** argv) +{ + int rank = 0, size = 0, buf = 0; + uint32_t ySize = 0, xSize = 0; + double* arr = 0; - calc(arr, ySize, xSize, rank, size); + MPI_Init(&argc, &argv); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); - if (rank == 0) - { - // Prepare output file - std::ofstream output(argv[2]); - if (!output.is_open()) + if (rank == 0) { - std::cout << "[Error] Can't open " << argv[2] << " for read\n"; - delete arr; - return 1; + // Check arguments + if (argc != 3) + { + std::cout << "[Error] Usage \n"; + buf = 1; + MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); + return 1; + } + + // Prepare input file + std::ifstream input(argv[1]); + if (!input.is_open()) + { + std::cout << "[Error] Can't open " << argv[1] << " for write\n"; + buf = 1; + MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); + return 1; + } + + // Read arguments from input + input >> ySize >> xSize; + MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); + + arr = new double[ySize * xSize]; + + for (uint32_t y = 0; y < ySize; y++) + { + for (uint32_t x = 0; x < xSize; x++) + { + input >> arr[y*xSize + x]; + } + } + input.close(); + } else { + MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); + if (buf != 0) + { + return 1; + } } - for (uint32_t y = 0; y < ySize; y++) + calc(arr, ySize, xSize, rank, size); + + if (rank == 0) { - for (uint32_t x = 0; x < xSize; x++) - { - output << " " << arr[y*xSize + x]; - } - output << std::endl; + // Prepare output file + std::ofstream output(argv[2]); + if (!output.is_open()) + { + std::cout << "[Error] Can't open " << argv[2] << " for read\n"; + delete arr; + return 1; + } + for (uint32_t y = 0; y < ySize; y++) + { + for (uint32_t x = 0; x < xSize; x++) + { + output << " " << arr[y*xSize + x]; + } + output << std::endl; + } + output.close(); + delete arr; } - output.close(); - delete arr; - } - MPI_Finalize(); - return 0; + MPI_Finalize(); + return 0; } diff --git a/Loop2/src/main.cpp b/Loop2/src/main.cpp index 51ea220..bf99524 100644 --- a/Loop2/src/main.cpp +++ b/Loop2/src/main.cpp @@ -5,18 +5,91 @@ #include #include -void calc(double* arr, uint32_t ySize, uint32_t xSize, int rank, int size) -{ - if (rank == 0 && size > 0) - { - for (uint32_t y = 0; y < ySize - 1; y++) - { - for (uint32_t x = 3; x < xSize; x++) - { - arr[y*xSize + x] = sin(0.00001*arr[(y + 1)*xSize + x - 3]); - } +#include + +#define ROOT 0 +// auntidep +// y = 1; y < ySize; +// x = 0; y < xSize - 3 +//arr[(y - 1)*xSize + (x + 3)] = sin(0.00001*arr[y*xSize + x]); + +void calc(double* arr, uint32_t ySize, uint32_t xSize, int rank, int size) { + int* send_dist = NULL; + int* send_range = NULL; + + int* recv_dist = NULL; + int* recv_range = NULL; + + MPI_Bcast(&ySize, 1, MPI_UNSIGNED, ROOT, MPI_COMM_WORLD); + MPI_Bcast(&xSize, 1, MPI_UNSIGNED, ROOT, MPI_COMM_WORLD); + int real_rank_size = 0; + + if(rank == ROOT) { + size_t my_calc_size = ceil((ySize) * 1.0 / size); + // if proc_rank > ceil(ySize / my_calc_size) it will do nothing + real_rank_size = ceil((ySize) * 1.0 / my_calc_size); + size_t now_size = 0; + + recv_dist = (int*)calloc(size, sizeof(int)); + recv_range = (int*)calloc(size, sizeof(int)); + + send_dist = (int*)calloc(size, sizeof(int)); + send_range = (int*)calloc(size, sizeof(int)); + + for(int i = 0; i < real_rank_size; ++i) { + recv_dist[i] = i * xSize * my_calc_size; + recv_range[i] = xSize * std::min(my_calc_size, ySize - now_size); + now_size += std::min(my_calc_size, ySize - now_size); + } + + for(int i = 0; i < real_rank_size; ++i) { + send_range[i] = (i == real_rank_size - 1)? recv_range[i] : recv_range[i] + xSize; + send_dist[i] = recv_dist[i]; + } + } + int send_step = 0; + int recv_step = 0; + + + // Bug with MPI <- i hate it! really! + MPI_Scatter( + send_range, 1, MPI_INT, + &send_step, size, MPI_INT, + ROOT, MPI_COMM_WORLD); + MPI_Scatter( + recv_range, 1, MPI_INT, + &recv_step, size, MPI_INT, + ROOT, MPI_COMM_WORLD); + + double* my_calc = NULL; + my_calc = (double*)calloc(send_step, sizeof(double)); + + MPI_Scatterv( + arr, send_range, send_dist, MPI_DOUBLE, + my_calc, send_step, MPI_DOUBLE, + ROOT, MPI_COMM_WORLD); + + if((uint32_t)send_step > xSize) { + for(size_t y = 0; y < send_step / xSize - 1; ++y) { + for(size_t x = 3; x < xSize; ++x) { + my_calc[y*xSize + x] = sin(0.00001*my_calc[(y + 1) * xSize + (x - 3)]); + } + } + } + + MPI_Gatherv( + my_calc, recv_step, MPI_DOUBLE, + arr, recv_range, recv_dist, MPI_DOUBLE, + ROOT, MPI_COMM_WORLD); + + + free(my_calc); + if(rank == ROOT) { + free(recv_dist); + free(recv_range); + free(send_dist); + free(send_range); } - } } int main(int argc, char** argv) diff --git a/Loop3/src/main.cpp b/Loop3/src/main.cpp index 43cce1f..0d9c88f 100644 --- a/Loop3/src/main.cpp +++ b/Loop3/src/main.cpp @@ -5,18 +5,76 @@ #include #include +#define ROOT 0 + void calc(double* arr, uint32_t ySize, uint32_t xSize, int rank, int size) { - if (rank == 0 && size > 0) - { - for (uint32_t y = 4; y < ySize; y++) - { - for (uint32_t x = 0; x < xSize; x++) - { - arr[y*xSize + x] = sin(arr[(y - 4)*xSize + x]); - } + double* new_arr = NULL; + if(rank == ROOT) { + new_arr = (double*)calloc(xSize * ySize, sizeof(double)); + for(size_t y = 0; y < ySize; ++y) { + for(size_t x = 0; x < xSize; ++x) { + new_arr[x*ySize + y] = arr[y*xSize + x]; + } + } } - } + + int32_t step = 0; + int32_t* range = NULL; + int32_t* displace = NULL; + + MPI_Bcast(&xSize, 1, MPI_UNSIGNED, ROOT, MPI_COMM_WORLD); + MPI_Bcast(&ySize, 1, MPI_UNSIGNED, ROOT, MPI_COMM_WORLD); + + if(rank == ROOT) { + uint32_t x_step = ceil(xSize / (1.0 * size)); + uint32_t prev_range = 0; + + range = (int32_t*)calloc(size, sizeof(int32_t)); + displace = (int32_t*)calloc(size, sizeof(int32_t)); + + for(int send_rank = 0; send_rank < size; ++send_rank) { + // We will send start_y, step, xSize and part of arr + uint32_t real_step = std::min(x_step, xSize - prev_range); + range[send_rank] = real_step * ySize; + displace[send_rank] = prev_range * ySize; + prev_range += real_step; + } + } + MPI_Scatter(range, 1, MPI_INT, + &step, 1, MPI_INT, + ROOT, MPI_COMM_WORLD); + + double* my_copy = NULL; + my_copy = (double*) calloc(step, sizeof(double)); + MPI_Scatterv(new_arr, range, displace, MPI_DOUBLE, + my_copy, step, MPI_DOUBLE, + ROOT, MPI_COMM_WORLD); + for(size_t x = 0; x < step / ySize; ++x) { + for(size_t y = 4; y < ySize; ++y) { + my_copy[x*ySize + y] = sin(my_copy[x*ySize + y - 4]); + } + } + + MPI_Gatherv(my_copy, step, MPI_DOUBLE, + new_arr, range, displace, MPI_DOUBLE, + ROOT, MPI_COMM_WORLD); + + if(rank == ROOT) { + for(size_t y = 0; y < ySize; ++y) { + for(size_t x = 0; x < xSize; ++x) { + arr[y*xSize + x] = new_arr[x*ySize + y]; + } + } + } + + + if(rank == ROOT) { + free(new_arr); + free(range); + free(displace); + } + free(my_copy); } int main(int argc, char** argv) diff --git a/Loop4/src/main.cpp b/Loop4/src/main.cpp index d958d54..71e7607 100644 --- a/Loop4/src/main.cpp +++ b/Loop4/src/main.cpp @@ -5,93 +5,177 @@ #include #include +#define ROOT 0 + +// arr[z*ySize*xSize + y*xSize + x] = sin(arr[(z - 1)*ySize*xSize + (y + 1)*xSize + x + 1]); void calc(double* arr, uint32_t zSize, uint32_t ySize, uint32_t xSize, int rank, int size) { - if (rank == 0 && size > 0) { - for (uint32_t z = 1; z < zSize; z++) { - for (uint32_t y = 0; y < ySize - 1; y++) { - for (uint32_t x = 0; x < xSize - 1; x++) { - arr[z*ySize*xSize + y*xSize + x] = sin(arr[(z - 1)*ySize*xSize + (y + 1)*xSize + x + 1]); + int* send_dist = NULL; + int* send_range = NULL; + + int* recv_dist = NULL; + int* recv_range = NULL; + + double* new_arr = NULL; + + MPI_Bcast(&zSize, 1, MPI_UNSIGNED, ROOT, MPI_COMM_WORLD); + MPI_Bcast(&ySize, 1, MPI_UNSIGNED, ROOT, MPI_COMM_WORLD); + MPI_Bcast(&xSize, 1, MPI_UNSIGNED, ROOT, MPI_COMM_WORLD); + int real_rank_size = 0; + + if(rank == ROOT) { + size_t my_calc_size = ceil((ySize) * 1.0 / size); + // if proc_rank > ceil(ySize / my_calc_size) it will do nothing + real_rank_size = ceil((ySize) * 1.0 / my_calc_size); + size_t now_size = 0; + + new_arr = (double*)calloc(xSize*ySize, sizeof(double)); + + recv_dist = (int*)calloc(size, sizeof(int)); + recv_range = (int*)calloc(size, sizeof(int)); + + send_dist = (int*)calloc(size, sizeof(int)); + send_range = (int*)calloc(size, sizeof(int)); + + for(int i = 0; i < real_rank_size; ++i) { + recv_dist[i] = i * xSize * my_calc_size; + recv_range[i] = xSize * std::min(my_calc_size, ySize - now_size); + now_size += std::min(my_calc_size, ySize - now_size); + } + + for(int i = 0; i < real_rank_size; ++i) { + send_range[i] = (i == real_rank_size - 1)? recv_range[i] : recv_range[i] + xSize; + send_dist[i] = recv_dist[i]; + } + } + int send_step = 0; + int recv_step = 0; + + + // Bug with MPI <- i hate it! really! + MPI_Scatter( + send_range, 1, MPI_INT, + &send_step, size, MPI_INT, + ROOT, MPI_COMM_WORLD); + MPI_Scatter( + recv_range, 1, MPI_INT, + &recv_step, size, MPI_INT, + ROOT, MPI_COMM_WORLD); + + double* my_calc = NULL; + my_calc = (double*)calloc(send_step, sizeof(double)); + + for(size_t z = 0; z < zSize - 1; ++z) { + MPI_Scatterv( + arr + z*xSize*ySize, send_range, send_dist, MPI_DOUBLE, + my_calc, send_step, MPI_DOUBLE, + ROOT, MPI_COMM_WORLD); + + if((uint32_t)send_step > xSize) { + for(size_t y = 0; y < send_step / xSize - 1; ++y) { + for(size_t x = 0; x < xSize - 1; ++x) { + my_calc[y*xSize + x] = sin(my_calc[(y + 1) * xSize + (x + 1)]); + } + } + } + + MPI_Gatherv( + my_calc, recv_step, MPI_DOUBLE, + new_arr, recv_range, recv_dist, MPI_DOUBLE, + ROOT, MPI_COMM_WORLD); + if(rank == ROOT) { + for(size_t y = 0; y < xSize - 1; ++y) { + for(size_t x = 0; x < xSize - 1; ++x) { + (arr + (z + 1)*xSize*ySize)[y * xSize + x] = new_arr[y * xSize + x]; + } + } } - } } - } + + free(my_calc); + if(rank == ROOT) { + free(new_arr); + free(recv_dist); + free(recv_range); + free(send_dist); + free(send_range); + } } int main(int argc, char** argv) { - int rank = 0, size = 0, buf = 0; - uint32_t zSize = 0, ySize = 0, xSize = 0; - double* arr = 0; - - MPI_Init(&argc, &argv); - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &size); - - if (rank == 0) - { - // Check arguments - if (argc != 3) - { - std::cout << "[Error] Usage \n"; - buf = 1; - MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); - return 1; - } + int rank = 0, size = 0, buf = 0; + uint32_t zSize = 0, ySize = 0, xSize = 0; + double* arr = 0; - // Prepare input file - std::ifstream input(argv[1]); - if (!input.is_open()) + MPI_Init(&argc, &argv); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + + if (rank == 0) { - std::cout << "[Error] Can't open " << argv[1] << " for write\n"; - buf = 1; - MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); - return 1; - } + // Check arguments + if (argc != 3) + { + std::cout << "[Error] Usage \n"; + buf = 1; + MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); + return 1; + } + + // Prepare input file + std::ifstream input(argv[1]); + if (!input.is_open()) + { + std::cout << "[Error] Can't open " << argv[1] << " for write\n"; + buf = 1; + MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); + return 1; + } - // Read arguments from input - input >> zSize >> ySize >> xSize; - MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); + // Read arguments from input + input >> zSize >> ySize >> xSize; + MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); - arr = new double[zSize * ySize * xSize]; - for (uint32_t z = 0; z < zSize; z++) { - for (uint32_t y = 0; y < ySize; y++) { - for (uint32_t x = 0; x < xSize; x++) { - input >> arr[z*ySize*xSize + y*xSize + x]; + arr = new double[zSize * ySize * xSize]; + for (uint32_t z = 0; z < zSize; z++) { + for (uint32_t y = 0; y < ySize; y++) { + for (uint32_t x = 0; x < xSize; x++) { + input >> arr[z*ySize*xSize + y*xSize + x]; + } + } + } + input.close(); + } else { + MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); + if (buf != 0) + { + return 1; } - } - } - input.close(); - } else { - MPI_Bcast(&buf, 1, MPI_INT, 0, MPI_COMM_WORLD); - if (buf != 0) - { - return 1; } - } - calc(arr, zSize, ySize, xSize, rank, size); + calc(arr, zSize, ySize, xSize, rank, size); - if (rank == 0) - { - // Prepare output file - std::ofstream output(argv[2]); - if (!output.is_open()) + if (rank == 0) { - std::cout << "[Error] Can't open " << argv[2] << " for read\n"; - delete arr; - return 1; - } + // Prepare output file + std::ofstream output(argv[2]); + if (!output.is_open()) + { + std::cout << "[Error] Can't open " << argv[2] << " for read\n"; + delete arr; + return 1; + } - for (uint32_t z = 0; z < zSize; z++) { - for (uint32_t y = 0; y < ySize; y++) { - for (uint32_t x = 0; x < xSize; x++) { - output << " " << arr[z*ySize*xSize + y*xSize + x]; + for (uint32_t z = 0; z < zSize; z++) { + for (uint32_t y = 0; y < ySize; y++) { + for (uint32_t x = 0; x < xSize; x++) { + output << " " << arr[z*ySize*xSize + y*xSize + x]; + } + output << std::endl; + } + output << std::endl; } - output << std::endl; - } - output << std::endl; - } output.close(); delete arr; } diff --git a/OpenMP2/src/main.cpp b/OpenMP2/src/main.cpp index 8dab215..cfe430c 100644 --- a/OpenMP2/src/main.cpp +++ b/OpenMP2/src/main.cpp @@ -3,9 +3,20 @@ #include #include -double calc(uint32_t x_last, uint32_t num_threads) -{ - return 0; +double calc(uint32_t x_last, uint32_t num_threads) { + double res = 0; + double* res_buf = (double*)calloc(x_last, sizeof(double)); + #pragma omp parallel num_threads(num_threads) + { + #pragma omp for + for(int i = x_last; i > 0; --i) { + res_buf[i - 1] = 1.0 / i; + } + } + for(int i = x_last; i > 0; --i) { + res += res_buf[i - 1]; + } + return res; } int main(int argc, char** argv) diff --git a/OpenMP3/src/main.cpp b/OpenMP3/src/main.cpp index 416ee06..4c20380 100644 --- a/OpenMP3/src/main.cpp +++ b/OpenMP3/src/main.cpp @@ -2,16 +2,43 @@ #include #include #include -#include +#include -double func(double x) -{ - return sin(x); -} +#define PI (2 * asin(1)) +#define DELTA(dx) (PI / dx) + +double (*f)(double) = sin; double calc(double x0, double x1, double dx, uint32_t num_threads) -{ - return 0; +{ + if(x1 - x0 == 0) { + return 0.0; + } + double res = 0.0; + size_t n = size_t((x1 - x0) / dx) + 1; // #pragma omp for don't work with double type + double* resbuf = (double*)calloc(n, sizeof(double)); + #pragma omp parallel num_threads(num_threads) + { + #pragma omp for + for(size_t _x = 0; _x < n; ++_x) { + if(_x == 0 || _x == n - 1) { + resbuf[_x] = f(x0 + _x * dx) / 2 * dx; + } else { + resbuf[_x ]= f(x0 + _x * dx) * dx; + } + } + } + size_t j = 0; + size_t s = 0; + for(size_t i = 0; i <= n; ++i) { + res += resbuf[j]; + j += DELTA(dx); + if(j >= n) { + s++; + j = s; + } + } + return res; } int main(int argc, char** argv) diff --git a/OpenMP4/run.sh b/OpenMP4/run.sh index 6b269de..39f9d58 100755 --- a/OpenMP4/run.sh +++ b/OpenMP4/run.sh @@ -1,7 +1,11 @@ #! /bin/bash # run.sh compiler="g++" +<<<<<<< HEAD +flags="-fopenmp -lm -std=c++11" +======= flags="-fopenmp -std=c++11 -Wall -Wextra -Werror" +>>>>>>> 70a4e91123320376c0d66f4a55d8449031b08bfd src="./src/main.cpp" build="./build" exe="$build/task" diff --git a/OpenMP4/src/main.cpp b/OpenMP4/src/main.cpp index 036d8e3..bc38cd3 100644 --- a/OpenMP4/src/main.cpp +++ b/OpenMP4/src/main.cpp @@ -2,10 +2,106 @@ #include #include #include +#include + +#include +/* +typedef struct { + double value; + size_t gr; +} limd_t; + +limd_t limd_add(limd_t lv, limd_t rv) { + size_t new_gr = 0; + if(lv.gr > rv.gr) { + new_gr = lv.gr = rv.gr; + lv.value /= pow(10, lv.gr - rv.gr); + } else { + new_gr = rv.gr = lv.gr; + lv.value /= pow(10, rv.gr - lv.gr); + } + limd_t res; + res.value = lv.value + rv.value; + res.gr = new_gr; + return res; +} + +limd_t limd_mul(limd_t lv, limd_t rv) { + limd_t res; + res.value = lv.value * rv.value; + res.gr = lv.gr + rv.gr; + return res; +} + +limd_t limd_sub(limd_t lv, limd_t rv) { + rv.value *= -1; + return limd_add(lv, rv); +} + +limd_t limd_dev(limd_t lv, limd_t rv) { + rv.value = 1 / rv.value; + rv.gr *= -1; + return limd_mul(lv, rv); +} + +limd_t limd_init(double value) { + if(value == 0) { + return {0, (size_t(-1))}; + } + limd_t res = {value, 0}; + while(res.value > 10 && res.value < 1) { + if(res.value > 1.0) { + res.value /= 10; + res.gr--; + } else { + res.value *= 10; + res.gr++; + } + } + return res; +}*/ double calc(uint32_t x_last, uint32_t num_threads) { - return 0; + double exp = 0.0; + std::vector fact(x_last); + for(auto i = next(fact.begin()); i != fact.end(); i = next(i)) { + *i = -1; + } + fact[0] = 1; + double c = 0.0; + #pragma omp parallel num_threads(num_threads) + { + #pragma omp for + for(size_t n = 1; n <= x_last; ++n) { + double my_fact = ( 1.0 / n ); + for(size_t i = n - 1; i >= 0; --i) { + if(fact[i] == -1) { + my_fact *= ( 1.0 / i ); + } else if(fact[i] > 0) { + my_fact *= fact[i]; + fact[n] = my_fact; + break; + } else { // fact[i] < 0 + break; + } + } + //std::cout << "I HAVE " << my_fact << std::endl; + /* + #pragma omp critical + { + double y = my_fact - c; + double t = exp + y; + c = (t - exp) - y; + exp = t; + }*/ + } + } + + for(auto i = fact.rbegin(); i != fact.rend(); i = next(i)) { + exp += *i; + } + return exp; } int main(int argc, char** argv)