diff --git a/geometry/graham_scan_algorithm.cpp b/geometry/graham_scan_algorithm.cpp index 06c1c0ec6..916f244d7 100644 --- a/geometry/graham_scan_algorithm.cpp +++ b/geometry/graham_scan_algorithm.cpp @@ -1,61 +1,58 @@ /****************************************************************************** -* @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 IO Operations -#include /// for std::assert -#include /// for std::vector -#include "./graham_scan_functions.hpp" /// for all the functions used + * @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::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 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; diff --git a/geometry/graham_scan_functions.hpp b/geometry/graham_scan_functions.hpp index 9d9325289..976d8ad4c 100644 --- a/geometry/graham_scan_functions.hpp +++ b/geometry/graham_scan_functions.hpp @@ -1,44 +1,48 @@ /****************************************************************************** -* @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 IO operations -#include /// for std::stack -#include /// for std::vector -#include /// for std::swap -#include /// for mathematics and datatype conversion + * @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::grahamscan @@ -46,148 +50,150 @@ *******************************************************************************/ namespace geometry::grahamscan { - /****************************************************************************** - * @struct Point - * @brief for X and Y co-ordinates of the co-ordinate. - *******************************************************************************/ - struct Point { - int x, y; - }; +/****************************************************************************** + * @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; +// 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 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; } - /****************************************************************************** - * @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); + 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; + } } - /****************************************************************************** - * @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); + // Place the bottom-most point at first position + std::swap(points[0], points[min]); - if (val == 0) return 0; // collinear - return (val > 0) ? 1 : 2; // clock or counter-clock wise + // 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 } - /****************************************************************************** - * @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); + // If modified array of points has less than 3 points, convex hull is not + // possible + if (m < 3) + return {}; - // Find orientation - int o = orientation(p0, *p1, *p2); - if (o == 0) { - return (distSq(p0, *p2) >= distSq(p0, *p1)) ? -1 : 1; - } + // Create an empty stack and push first three points to it. + std::stack S; + S.push(points[0]); + S.push(points[1]); + S.push(points[2]); - 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 S; - S.push(points[0]); - S.push(points[1]); - S.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 (S.size() > 1 && orientation(nextToTop(S), S.top(), points[i]) != 2) { - S.pop(); - } - S.push(points[i]); - } - - std::vector result; - // Now stack has the output points, push them into the resultant vector - while (!S.empty()) { - Point p = S.top(); - result.push_back(p); + // 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 (S.size() > 1 && + orientation(nextToTop(S), S.top(), points[i]) != 2) { S.pop(); } - - return result; // return resultant vector with Convex Hull co-ordinates. + S.push(points[i]); } + + std::vector result; + // Now stack has the output points, push them into the resultant vector + while (!S.empty()) { + Point p = S.top(); + result.push_back(p); + S.pop(); + } + + return result; // return resultant vector with Convex Hull co-ordinates. +} } // namespace geometry::grahamscan diff --git a/numerical_methods/composite_simpson_rule.cpp b/numerical_methods/composite_simpson_rule.cpp index 9a5e8c180..085b431ff 100644 --- a/numerical_methods/composite_simpson_rule.cpp +++ b/numerical_methods/composite_simpson_rule.cpp @@ -43,6 +43,8 @@ #include /// for IO operations #include /// for std::map container +#include "math.h" + /** * @namespace numerical_methods * @brief Numerical algorithms/methods @@ -64,13 +66,13 @@ namespace simpson_method { * @returns the result of the integration */ double evaluate_by_simpson(std::int32_t N, double h, double a, - std::function func) { + const std::function& func) { std::map data_table; // Contains the data points. key: i, value: f(xi) double xi = a; // Initialize xi to the starting point x0 = a // Create the data table - double temp; + double temp = NAN; for (std::int32_t i = 0; i <= N; i++) { temp = func(xi); data_table.insert( @@ -82,12 +84,13 @@ double evaluate_by_simpson(std::int32_t N, double h, double a, // Remember: f(x0) + 4*f(x1) + 2*f(x2) + ... + 2*f(xN-2) + 4*f(xN-1) + f(xN) double evaluate_integral = 0; for (std::int32_t i = 0; i <= N; i++) { - if (i == 0 || i == N) + if (i == 0 || i == N) { evaluate_integral += data_table.at(i); - else if (i % 2 == 1) + } else if (i % 2 == 1) { evaluate_integral += 4 * data_table.at(i); - else + } else { evaluate_integral += 2 * data_table.at(i); + } } // Multiply by the coefficient h/3 @@ -170,7 +173,7 @@ int main(int argc, char** argv) { /// interval. MUST BE EVEN double a = 1, b = 3; /// Starting and ending point of the integration in /// the real axis - double h; /// Step, calculated by a, b and N + double h = NAN; /// Step, calculated by a, b and N bool used_argv_parameters = false; // If argv parameters are used then the assert must be omitted @@ -180,18 +183,20 @@ int main(int argc, char** argv) { // displaying messages) if (argc == 4) { N = std::atoi(argv[1]); - a = (double)std::atof(argv[2]); - b = (double)std::atof(argv[3]); + a = std::atof(argv[2]); + b = std::atof(argv[3]); // Check if a 0 && "N has to be > 0"); - if (N < 16 || a != 1 || b != 3) + if (N < 16 || a != 1 || b != 3) { used_argv_parameters = true; + } std::cout << "You selected N=" << N << ", a=" << a << ", b=" << b << std::endl; - } else + } else { std::cout << "Default N=" << N << ", a=" << a << ", b=" << b << std::endl; + } // Find the step h = (b - a) / N; diff --git a/numerical_methods/inverse_fast_fourier_transform.cpp b/numerical_methods/inverse_fast_fourier_transform.cpp index d2248be7b..0970d40cd 100644 --- a/numerical_methods/inverse_fast_fourier_transform.cpp +++ b/numerical_methods/inverse_fast_fourier_transform.cpp @@ -4,10 +4,10 @@ * (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). + * 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) */ @@ -23,14 +23,15 @@ */ namespace numerical_methods { /** - * @brief InverseFastFourierTransform is a recursive function which returns list of - * complex numbers + * @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) { +std::complex *InverseFastFourierTransform(std::complex *p, + uint8_t n) { if (n == 1) { return p; /// Base Case To return } @@ -39,9 +40,9 @@ std::complex *InverseFastFourierTransform(std::complex *p, uint8 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 + + 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 @@ -52,8 +53,9 @@ std::complex *InverseFastFourierTransform(std::complex *p, uint8 if (j % 2 == 0) { pe[k1++] = p[j]; /// Assigning values of even Coefficients - } else + } else { po[k2++] = p[j]; /// Assigning value of odd Coefficients + } } std::complex *ye = @@ -75,12 +77,10 @@ std::complex *InverseFastFourierTransform(std::complex *p, uint8 k1++; k2++; } - - if(n!=2){ - + + if (n != 2) { delete[] pe; delete[] po; - } delete[] ye; /// Deleting dynamic array ye @@ -118,16 +118,17 @@ static void test() { 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); + 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++) { @@ -135,10 +136,8 @@ static void test() { (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;