From ea76786f1274be7e2686dc7ecd45d9850b3203f1 Mon Sep 17 00:00:00 2001 From: Lajat5 <64376519+Lazeeez@users.noreply.github.com> Date: Sun, 14 Nov 2021 22:26:46 +0530 Subject: [PATCH] feat: Added Graham Scan Algorithm. (#1836) * Implementated Grahamscan Algorithm for Convex Hull * Update graham_scan_algorithm.cpp * Update graham_scan_functions.h * Update graham_scan_algorithm.cpp * Update graham_scan_functions.h * Update graham_scan_algorithm.cpp * Update and rename graham_scan_functions.h to graham_scan_functions.hpp * Update geometry/graham_scan_algorithm.cpp Co-authored-by: David Leal * Update geometry/graham_scan_algorithm.cpp Co-authored-by: David Leal * Update geometry/graham_scan_functions.hpp Co-authored-by: David Leal * Update geometry/graham_scan_functions.hpp Co-authored-by: David Leal * updating DIRECTORY.md * clang-format and clang-tidy fixes for e89e4c8c * clang-format and clang-tidy fixes for 7df4778f * Fix #1 * Update graham_scan_functions.hpp * Delete composite_simpson_rule.cpp * Delete inverse_fast_fourier_transform.cpp * Fix #2 * updating DIRECTORY.md * clang-format and clang-tidy fixes for 69b6832b * Fix #3 * clang-format and clang-tidy fixes for 1c05ca7c * Update graham_scan_functions.hpp * Fix #4 * clang-format and clang-tidy fixes for 2957fd21 * Create composite_simpson_rule.cpp * updating DIRECTORY.md * Create inverse_fast_fourier_transform.cpp * updating DIRECTORY.md * clang-format and clang-tidy fixes for 405d21a5 * clang-format and clang-tidy fixes for 333ef5ca * Update geometry/graham_scan_functions.hpp Co-authored-by: David Leal * Update geometry/graham_scan_algorithm.cpp Co-authored-by: David Leal * Update geometry/graham_scan_algorithm.cpp Co-authored-by: David Leal * Update geometry/graham_scan_algorithm.cpp Co-authored-by: David Leal * Update geometry/graham_scan_functions.hpp Co-authored-by: David Leal * Update geometry/graham_scan_algorithm.cpp Co-authored-by: David Leal * Update geometry/graham_scan_algorithm.cpp Co-authored-by: David Leal * clang-format and clang-tidy fixes for ee4cb635 * Update graham_scan_algorithm.cpp * Update graham_scan_functions.hpp * clang-format and clang-tidy fixes for f2f69234 * Update graham_scan_functions.hpp * Create partition_problem.cpp * Update partition_problem.cpp * Delete partition_problem.cpp Co-authored-by: David Leal Co-authored-by: github-actions <${GITHUB_ACTOR}@users.noreply.github.com> --- DIRECTORY.md | 2 + geometry/graham_scan_algorithm.cpp | 76 +++++ geometry/graham_scan_functions.hpp | 209 ++++++++++++ .../inverse_fast_fourier_transform.cpp | 320 +++++++++--------- 4 files changed, 447 insertions(+), 160 deletions(-) create mode 100644 geometry/graham_scan_algorithm.cpp create mode 100644 geometry/graham_scan_functions.hpp diff --git a/DIRECTORY.md b/DIRECTORY.md index 3ed3dfd60..392415722 100644 --- a/DIRECTORY.md +++ b/DIRECTORY.md @@ -109,6 +109,8 @@ * [Word Break](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/dynamic_programming/word_break.cpp) ## Geometry + * [Graham Scan Algorithm](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/geometry/graham_scan_algorithm.cpp) + * [Graham Scan Functions](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/geometry/graham_scan_functions.hpp) * [Jarvis Algorithm](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/geometry/jarvis_algorithm.cpp) * [Line Segment Intersection](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/geometry/line_segment_intersection.cpp) diff --git a/geometry/graham_scan_algorithm.cpp b/geometry/graham_scan_algorithm.cpp new file mode 100644 index 000000000..ac1435cef --- /dev/null +++ b/geometry/graham_scan_algorithm.cpp @@ -0,0 +1,76 @@ +/****************************************************************************** + * @file + * @brief Implementation of the [Convex + * Hull](https://en.wikipedia.org/wiki/Convex_hull) implementation using [Graham + * Scan](https://en.wikipedia.org/wiki/Graham_scan) + * @details + * In geometry, the convex hull or convex envelope or convex closure of a shape + * is the smallest convex set that contains it. The convex hull may be defined + * either as the intersection of all convex sets containing a given subset of a + * Euclidean space, or equivalently as the set of all convex combinations of + * points in the subset. For a bounded subset of the plane, the convex hull may + * be visualized as the shape enclosed by a rubber band stretched around the + * subset. + * + * The worst case time complexity of Jarvis’s Algorithm is O(n^2). Using + * Graham’s scan algorithm, we can find Convex Hull in O(nLogn) time. + * + * ### Implementation + * + * Sort points + * We first find the bottom-most point. The idea is to pre-process + * points be sorting them with respect to the bottom-most point. Once the points + * are sorted, they form a simple closed path. + * The sorting criteria is to use the orientation to compare angles without + * actually computing them (See the compare() function below) because + * computation of actual angles would be inefficient since trigonometric + * functions are not simple to evaluate. + * + * Accept or Reject Points + * Once we have the closed path, the next step is to traverse the path and + * remove concave points on this path using orientation. The first two points in + * sorted array are always part of Convex Hull. For remaining points, we keep + * track of recent three points, and find the angle formed by them. Let the + * three points be prev(p), curr(c) and next(n). If the orientation of these + * points (considering them in the same order) is not counterclockwise, we + * discard c, otherwise we keep it. + * + * @author [Lajat Manekar](https://github.com/Lazeeez) + * + *******************************************************************************/ +#include /// for std::assert +#include /// for IO Operations +#include /// for std::vector + +#include "./graham_scan_functions.hpp" /// for all the functions used + +/******************************************************************************* + * @brief Self-test implementations + * @returns void + *******************************************************************************/ +static void test() { + std::vector points = { + {0, 3}, {1, 1}, {2, 2}, {4, 4}, {0, 0}, {1, 2}, {3, 1}, {3, 3}}; + std::vector expected_result = { + {0, 3}, {4, 4}, {3, 1}, {0, 0}}; + std::vector derived_result; + std::vector res; + + derived_result = geometry::grahamscan::convexHull(points, points.size()); + + std::cout << "1st test: "; + for (int i = 0; i < expected_result.size(); i++) { + assert(derived_result[i].x == expected_result[i].x); + assert(derived_result[i].y == expected_result[i].y); + } + std::cout << "passed!" << std::endl; +} + +/******************************************************************************* + * @brief Main function + * @returns 0 on exit + *******************************************************************************/ +int main() { + test(); // run self-test implementations + return 0; +} diff --git a/geometry/graham_scan_functions.hpp b/geometry/graham_scan_functions.hpp new file mode 100644 index 000000000..f6e05095e --- /dev/null +++ b/geometry/graham_scan_functions.hpp @@ -0,0 +1,209 @@ +/****************************************************************************** + * @file + * @brief Implementation of the [Convex + * Hull](https://en.wikipedia.org/wiki/Convex_hull) implementation using [Graham + * Scan](https://en.wikipedia.org/wiki/Graham_scan) + * @details + * In geometry, the convex hull or convex envelope or convex closure of a shape + * is the smallest convex set that contains it. The convex hull may be defined + * either as the intersection of all convex sets containing a given subset of a + * Euclidean space, or equivalently as the set of all convex combinations of + * points in the subset. For a bounded subset of the plane, the convex hull may + * be visualized as the shape enclosed by a rubber band stretched around the + * subset. + * + * The worst case time complexity of Jarvis’s Algorithm is O(n^2). Using + * Graham’s scan algorithm, we can find Convex Hull in O(nLogn) time. + * + * ### Implementation + * + * Sort points + * We first find the bottom-most point. The idea is to pre-process + * points be sorting them with respect to the bottom-most point. Once the points + * are sorted, they form a simple closed path. + * The sorting criteria is to use the orientation to compare angles without + * actually computing them (See the compare() function below) because + * computation of actual angles would be inefficient since trigonometric + * functions are not simple to evaluate. + * + * Accept or Reject Points + * Once we have the closed path, the next step is to traverse the path and + * remove concave points on this path using orientation. The first two points in + * sorted array are always part of Convex Hull. For remaining points, we keep + * track of recent three points, and find the angle formed by them. Let the + * three points be prev(p), curr(c) and next(n). If orientation of these points + * (considering them in same order) is not counterclockwise, we discard c, + * otherwise we keep it. + * + * @author [Lajat Manekar](https://github.com/Lazeeez) + * + *******************************************************************************/ +#include /// for std::swap +#include /// for mathematics and datatype conversion +#include /// for IO operations +#include /// for std::stack +#include /// for std::vector + +/****************************************************************************** + * @namespace geometry + * @brief geometric algorithms + *******************************************************************************/ +namespace geometry { + +/****************************************************************************** + * @namespace graham scan + * @brief convex hull algorithm + *******************************************************************************/ +namespace grahamscan { + +/****************************************************************************** + * @struct Point + * @brief for X and Y co-ordinates of the co-ordinate. + *******************************************************************************/ +struct Point { + int x, y; +}; + +// A global point needed for sorting points with reference +// to the first point Used in compare function of qsort() + +Point p0; + +/****************************************************************************** + * @brief A utility function to find next to top in a stack. + * @param S Stack to be used for the process. + * @returns @param Point Co-ordinates of the Point + *******************************************************************************/ +Point nextToTop(std::stack *S) { + Point p = S->top(); + S->pop(); + Point res = S->top(); + S->push(p); + return res; +} + +/****************************************************************************** + * @brief A utility function to return square of distance between p1 and p2. + * @param p1 Co-ordinates of Point 1 . + * @param p2 Co-ordinates of Point 2 . + * @returns @param int distance between p1 and p2. + *******************************************************************************/ +int distSq(Point p1, Point p2) { + return (p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y); +} + +/****************************************************************************** + * @brief To find orientation of ordered triplet (p, q, r). + * @param p Co-ordinates of Point p . + * @param q Co-ordinates of Point q . + * @param r Co-ordinates of Point r . + * @returns @param int 0 --> p, q and r are collinear, 1 --> Clockwise, + * 2 --> Counterclockwise + *******************************************************************************/ +int orientation(Point p, Point q, Point r) { + int val = (q.y - p.y) * (r.x - q.x) - (q.x - p.x) * (r.y - q.y); + + if (val == 0) { + return 0; // collinear + } + return (val > 0) ? 1 : 2; // clock or counter-clock wise +} + +/****************************************************************************** + * @brief A function used by library function qsort() to sort an array of + * points with respect to the first point + * @param vp1 Co-ordinates of Point 1 . + * @param vp2 Co-ordinates of Point 2 . + * @returns @param int distance between p1 and p2. + *******************************************************************************/ +int compare(const void *vp1, const void *vp2) { + auto *p1 = static_cast(vp1); + auto *p2 = static_cast(vp2); + + // Find orientation + int o = orientation(p0, *p1, *p2); + if (o == 0) { + return (distSq(p0, *p2) >= distSq(p0, *p1)) ? -1 : 1; + } + + return (o == 2) ? -1 : 1; +} + +/****************************************************************************** + * @brief Prints convex hull of a set of n points. + * @param points vector of Point with co-ordinates. + * @param size Size of the vector. + * @returns @param vector vector of Conver Hull. + *******************************************************************************/ +std::vector convexHull(std::vector points, uint64_t size) { + // Find the bottom-most point + int ymin = points[0].y, min = 0; + for (int i = 1; i < size; i++) { + int y = points[i].y; + + // Pick the bottom-most or chose the left-most point in case of tie + if ((y < ymin) || (ymin == y && points[i].x < points[min].x)) { + ymin = points[i].y, min = i; + } + } + + // Place the bottom-most point at first position + std::swap(points[0], points[min]); + + // Sort n-1 points with respect to the first point. A point p1 comes + // before p2 in sorted output if p2 has larger polar angle + // (in counterclockwise direction) than p1. + p0 = points[0]; + qsort(&points[1], size - 1, sizeof(Point), compare); + + // If two or more points make same angle with p0, Remove all but the one + // that is farthest from p0 Remember that, in above sorting, our criteria + // was to keep the farthest point at the end when more than one points have + // same angle. + int m = 1; // Initialize size of modified array + for (int i = 1; i < size; i++) { + // Keep removing i while angle of i and i+1 is same with respect to p0 + while (i < size - 1 && orientation(p0, points[i], points[i + 1]) == 0) { + i++; + } + + points[m] = points[i]; + m++; // Update size of modified array + } + + // If modified array of points has less than 3 points, convex hull is not + // possible + if (m < 3) { + return {}; + }; + + // Create an empty stack and push first three points to it. + std::stack St; + St.push(points[0]); + St.push(points[1]); + St.push(points[2]); + + // Process remaining n-3 points + for (int i = 3; i < m; i++) { + // Keep removing top while the angle formed by + // points next-to-top, top, and points[i] makes + // a non-left turn + while (St.size() > 1 && + orientation(nextToTop(&St), St.top(), points[i]) != 2) { + St.pop(); + } + St.push(points[i]); + } + + std::vector result; + // Now stack has the output points, push them into the resultant vector + while (!St.empty()) { + Point p = St.top(); + result.push_back(p); + St.pop(); + } + + return result; // return resultant vector with Convex Hull co-ordinates. +} +} // namespace grahamscan +} // namespace geometry diff --git a/numerical_methods/inverse_fast_fourier_transform.cpp b/numerical_methods/inverse_fast_fourier_transform.cpp index 0970d40cd..3837c21f9 100644 --- a/numerical_methods/inverse_fast_fourier_transform.cpp +++ b/numerical_methods/inverse_fast_fourier_transform.cpp @@ -1,160 +1,160 @@ -/** - * @file - * @brief [An inverse fast Fourier transform - * (IFFT)](https://www.geeksforgeeks.org/python-inverse-fast-fourier-transformation/) - * is an algorithm that computes the inverse fourier transform. - * @details - * This algorithm has an application in use case scenario where a user wants - * find coefficients of a function in a short time by just using points - * generated by DFT. Time complexity this algorithm computes the IDFT in - * O(nlogn) time in comparison to traditional O(n^2). - * @author [Ameya Chawla](https://github.com/ameyachawlaggsipu) - */ - -#include /// for assert -#include /// for mathematical-related functions -#include /// for storing points and coefficents -#include /// for IO operations -#include /// for std::vector - -/** - * @namespace numerical_methods - * @brief Numerical algorithms/methods - */ -namespace numerical_methods { -/** - * @brief InverseFastFourierTransform is a recursive function which returns list - * of complex numbers - * @param p List of Coefficents in form of complex numbers - * @param n Count of elements in list p - * @returns p if n==1 - * @returns y if n!=1 - */ -std::complex *InverseFastFourierTransform(std::complex *p, - uint8_t n) { - if (n == 1) { - return p; /// Base Case To return - } - - double pi = 2 * asin(1.0); /// Declaring value of pi - - std::complex om = std::complex( - cos(2 * pi / n), sin(2 * pi / n)); /// Calculating value of omega - - om.real(om.real() / n); /// One change in comparison with DFT - om.imag(om.imag() / n); /// One change in comparison with DFT - - auto *pe = new std::complex[n / 2]; /// Coefficients of even power - - auto *po = new std::complex[n / 2]; /// Coefficients of odd power - - int k1 = 0, k2 = 0; - for (int j = 0; j < n; j++) { - if (j % 2 == 0) { - pe[k1++] = p[j]; /// Assigning values of even Coefficients - - } else { - po[k2++] = p[j]; /// Assigning value of odd Coefficients - } - } - - std::complex *ye = - InverseFastFourierTransform(pe, n / 2); /// Recursive Call - - std::complex *yo = - InverseFastFourierTransform(po, n / 2); /// Recursive Call - - auto *y = new std::complex[n]; /// Final value representation list - - k1 = 0, k2 = 0; - - for (int i = 0; i < n / 2; i++) { - y[i] = - ye[k1] + pow(om, i) * yo[k2]; /// Updating the first n/2 elements - y[i + n / 2] = - ye[k1] - pow(om, i) * yo[k2]; /// Updating the last n/2 elements - - k1++; - k2++; - } - - if (n != 2) { - delete[] pe; - delete[] po; - } - - delete[] ye; /// Deleting dynamic array ye - delete[] yo; /// Deleting dynamic array yo - return y; -} - -} // namespace numerical_methods - -/** - * @brief Self-test implementations - * @details - * Declaring two test cases and checking for the error - * in predicted and true value is less than 0.000000000001. - * @returns void - */ -static void test() { - /* descriptions of the following test */ - - auto *t1 = new std::complex[2]; /// Test case 1 - auto *t2 = new std::complex[4]; /// Test case 2 - - t1[0] = {3, 0}; - t1[1] = {-1, 0}; - t2[0] = {10, 0}; - t2[1] = {-2, -2}; - t2[2] = {-2, 0}; - t2[3] = {-2, 2}; - - uint8_t n1 = 2; - uint8_t n2 = 4; - std::vector> r1 = { - {1, 0}, {2, 0}}; /// True Answer for test case 1 - - std::vector> r2 = { - {1, 0}, {2, 0}, {3, 0}, {4, 0}}; /// True Answer for test case 2 - - std::complex *o1 = - numerical_methods::InverseFastFourierTransform(t1, n1); - - std::complex *o2 = - numerical_methods::InverseFastFourierTransform(t2, n2); - - for (uint8_t i = 0; i < n1; i++) { - assert((r1[i].real() - o1[i].real() < 0.000000000001) && - (r1[i].imag() - o1[i].imag() < - 0.000000000001)); /// Comparing for both real and imaginary - /// values for test case 1 - } - - for (uint8_t i = 0; i < n2; i++) { - assert((r2[i].real() - o2[i].real() < 0.000000000001) && - (r2[i].imag() - o2[i].imag() < - 0.000000000001)); /// Comparing for both real and imaginary - /// values for test case 2 - } - - delete[] t1; - delete[] t2; - delete[] o1; - delete[] o2; - std::cout << "All tests have successfully passed!\n"; -} - -/** - * @brief Main function - * @param argc commandline argument count (ignored) - * @param argv commandline array of arguments (ignored) - * calls automated test function to test the working of fast fourier transform. - * @returns 0 on exit - */ - -int main(int argc, char const *argv[]) { - test(); // run self-test implementations - // with 2 defined test cases - return 0; -} +/** + * @file + * @brief [An inverse fast Fourier transform + * (IFFT)](https://www.geeksforgeeks.org/python-inverse-fast-fourier-transformation/) + * is an algorithm that computes the inverse fourier transform. + * @details + * This algorithm has an application in use case scenario where a user wants + * find coefficients of a function in a short time by just using points + * generated by DFT. Time complexity this algorithm computes the IDFT in + * O(nlogn) time in comparison to traditional O(n^2). + * @author [Ameya Chawla](https://github.com/ameyachawlaggsipu) + */ + +#include /// for assert +#include /// for mathematical-related functions +#include /// for storing points and coefficents +#include /// for IO operations +#include /// for std::vector + +/** + * @namespace numerical_methods + * @brief Numerical algorithms/methods + */ +namespace numerical_methods { +/** + * @brief InverseFastFourierTransform is a recursive function which returns list + * of complex numbers + * @param p List of Coefficents in form of complex numbers + * @param n Count of elements in list p + * @returns p if n==1 + * @returns y if n!=1 + */ +std::complex *InverseFastFourierTransform(std::complex *p, + uint8_t n) { + if (n == 1) { + return p; /// Base Case To return + } + + double pi = 2 * asin(1.0); /// Declaring value of pi + + std::complex om = std::complex( + cos(2 * pi / n), sin(2 * pi / n)); /// Calculating value of omega + + om.real(om.real() / n); /// One change in comparison with DFT + om.imag(om.imag() / n); /// One change in comparison with DFT + + auto *pe = new std::complex[n / 2]; /// Coefficients of even power + + auto *po = new std::complex[n / 2]; /// Coefficients of odd power + + int k1 = 0, k2 = 0; + for (int j = 0; j < n; j++) { + if (j % 2 == 0) { + pe[k1++] = p[j]; /// Assigning values of even Coefficients + + } else { + po[k2++] = p[j]; /// Assigning value of odd Coefficients + } + } + + std::complex *ye = + InverseFastFourierTransform(pe, n / 2); /// Recursive Call + + std::complex *yo = + InverseFastFourierTransform(po, n / 2); /// Recursive Call + + auto *y = new std::complex[n]; /// Final value representation list + + k1 = 0, k2 = 0; + + for (int i = 0; i < n / 2; i++) { + y[i] = + ye[k1] + pow(om, i) * yo[k2]; /// Updating the first n/2 elements + y[i + n / 2] = + ye[k1] - pow(om, i) * yo[k2]; /// Updating the last n/2 elements + + k1++; + k2++; + } + + if (n != 2) { + delete[] pe; + delete[] po; + } + + delete[] ye; /// Deleting dynamic array ye + delete[] yo; /// Deleting dynamic array yo + return y; +} + +} // namespace numerical_methods + +/** + * @brief Self-test implementations + * @details + * Declaring two test cases and checking for the error + * in predicted and true value is less than 0.000000000001. + * @returns void + */ +static void test() { + /* descriptions of the following test */ + + auto *t1 = new std::complex[2]; /// Test case 1 + auto *t2 = new std::complex[4]; /// Test case 2 + + t1[0] = {3, 0}; + t1[1] = {-1, 0}; + t2[0] = {10, 0}; + t2[1] = {-2, -2}; + t2[2] = {-2, 0}; + t2[3] = {-2, 2}; + + uint8_t n1 = 2; + uint8_t n2 = 4; + std::vector> r1 = { + {1, 0}, {2, 0}}; /// True Answer for test case 1 + + std::vector> r2 = { + {1, 0}, {2, 0}, {3, 0}, {4, 0}}; /// True Answer for test case 2 + + std::complex *o1 = + numerical_methods::InverseFastFourierTransform(t1, n1); + + std::complex *o2 = + numerical_methods::InverseFastFourierTransform(t2, n2); + + for (uint8_t i = 0; i < n1; i++) { + assert((r1[i].real() - o1[i].real() < 0.000000000001) && + (r1[i].imag() - o1[i].imag() < + 0.000000000001)); /// Comparing for both real and imaginary + /// values for test case 1 + } + + for (uint8_t i = 0; i < n2; i++) { + assert((r2[i].real() - o2[i].real() < 0.000000000001) && + (r2[i].imag() - o2[i].imag() < + 0.000000000001)); /// Comparing for both real and imaginary + /// values for test case 2 + } + + delete[] t1; + delete[] t2; + delete[] o1; + delete[] o2; + std::cout << "All tests have successfully passed!\n"; +} + +/** + * @brief Main function + * @param argc commandline argument count (ignored) + * @param argv commandline array of arguments (ignored) + * calls automated test function to test the working of fast fourier transform. + * @returns 0 on exit + */ + +int main(int argc, char const *argv[]) { + test(); // run self-test implementations + // with 2 defined test cases + return 0; +}