diff --git a/ComputationalGeometry/6Ex/PythonVersion/main.py b/ComputationalGeometry/6Ex/PythonVersion/main.py new file mode 100644 index 0000000..e1cee6b --- /dev/null +++ b/ComputationalGeometry/6Ex/PythonVersion/main.py @@ -0,0 +1,175 @@ +# Python3 program to find the minimum enclosing +# circle for N integer points in a 2-D plane +from math import sqrt +from random import randint, shuffle + +# Defining infinity +INF = 1e18 +MAX_POINT_COORD = 100 + + +# Structure to represent a 2D point +class Point: + def __init__(self, X=0, Y=0) -> None: + self.X = X + self.Y = Y + + +# Structure to represent a 2D circle +class Circle: + def __init__(self, c=Point(), r=0) -> None: + self.C = c + self.R = r + + def print(self) -> None: + print(f"(x - {self.C.X})^2 + (y - {self.C.Y})^2 = {self.R}^2") + + +# Function to return the euclidean distance +# between two points +def dist(a, b): + return sqrt(pow(a.X - b.X, 2) + + pow(a.Y - b.Y, 2)) + + +# Function to check whether a point lies inside +# or on the boundaries of the circle +def is_inside(c, p): + return dist(c.C, p) <= c.R + + +# The following two functions are used +# To find the equation of the circle when +# three points are given. + +# Helper method to get a circle defined by 3 points +def get_circle_center(bx, by, + cx, cy): + B = bx * bx + by * by + C = cx * cx + cy * cy + D = bx * cy - by * cx + return Point((cy * B - by * C) / (2 * D), + (bx * C - cx * B) / (2 * D)) + + +# Function to return the smallest circle +# that intersects 2 points +def circle_from1(A, B): + # Set the center to be the midpoint of A and B + C = Point((A.X + B.X) / 2.0, (A.Y + B.Y) / 2.0) + + # Set the radius to be half the distance AB + return Circle(C, dist(A, B) / 2.0) + + +# Function to return a unique circle that +# intersects three points +def circle_from2(A, B, C): + I = get_circle_center(B.X - A.X, B.Y - A.Y, + C.X - A.X, C.Y - A.Y) + + I.X += A.X + I.Y += A.Y + return Circle(I, dist(I, A)) + + +# Function to check whether a circle +# encloses the given points +def is_valid_circle(c, P): + # Iterating through all the points + # to check whether the points + # lie inside the circle or not + for p in P: + if (not is_inside(c, p)): + return False + return True + + +# Function to return the minimum enclosing +# circle for N <= 3 +def min_circle_trivial(P): + assert (len(P) <= 3) + if not P: + return Circle() + + elif (len(P) == 1): + return Circle(P[0], 0) + + elif (len(P) == 2): + return circle_from1(P[0], P[1]) + + # To check if MEC can be determined + # by 2 points only + for i in range(3): + for j in range(i + 1, 3): + + c = circle_from1(P[i], P[j]) + if (is_valid_circle(c, P)): + return c + + return circle_from2(P[0], P[1], P[2]) + + +# Returns the MEC using Welzl's algorithm +# Takes a set of input points P and a set R +# points on the circle boundary. +# n represents the number of points in P +# that are not yet processed. +def welzl_helper(P, R, n): + # Base case when all points processed or |R| = 3 + if (n == 0 or len(R) == 3): + return min_circle_trivial(R) + + # Pick a random point randomly + idx = randint(0, n - 1) + p = P[idx] + + # Put the picked point at the end of P + # since it's more efficient than + # deleting from the middle of the vector + P[idx], P[n - 1] = P[n - 1], P[idx] + + # Get the MEC circle d from the + # set of points P - :p + d = welzl_helper(P, R.copy(), n - 1) + + # If d contains p, return d + if (is_inside(d, p)): + return d + + # Otherwise, must be on the boundary of the MEC + R.append(p) + + # Return the MEC for P - :p and R U :p + return welzl_helper(P, R.copy(), n - 1) + + +def welzl(P): + P_copy = P.copy() + shuffle(P_copy) + return welzl_helper(P_copy, [], len(P_copy)) + + +def generate(num: int) -> list[Point]: + data: list[Point] = list() + + for i in range(num): + data.append(Point(randint(-MAX_POINT_COORD, MAX_POINT_COORD), randint(-MAX_POINT_COORD, MAX_POINT_COORD))) + + return data + + +def printPoints(data: list[Point]) -> None: + for i in range(len(data)): + print(f"{chr(i + 65)} = ({data[i].X}, {data[i].Y})") + + +# Driver code +if __name__ == '__main__': + num = int(input("Enter number of points: ")) + data = generate(num) + + printPoints(data) + + mec = welzl(data) + mec.print() diff --git a/ComputationalGeometry/6Ex/SCP.c b/ComputationalGeometry/6Ex/SCP.c index 1644590..1c2e749 100644 --- a/ComputationalGeometry/6Ex/SCP.c +++ b/ComputationalGeometry/6Ex/SCP.c @@ -24,21 +24,23 @@ double powd(double number) { circle primitive(point * N, int nlen) { if (nlen == 3) { - point temp; - circle crcl; - - for (int i = 0; i < nlen; ++i) { - temp = N[2]; - N[2] = N[i]; - N[i] = temp; - - crcl = centermass(N[0], N[1]); - if (belongs(crcl, N[2])) return crcl; - - N[i] = N[2]; - N[2] = temp; - } + for (int i = 0; i < 3; ++i) { + for (int j = i + 1; j < 3; ++j) { + circle crcl = centermass(N[i], N[j]); + bool all_inside = true; + for (int k = 0; k < 3; ++k) { + if (!belongs(crcl, N[k])) { + all_inside = false; + break; + } + } + if (all_inside) { + return crcl; + } + } + } + return byThreePoints(N); } else if (nlen == 2) { return centermass(N[0], N[1]); @@ -54,28 +56,28 @@ circle centermass(point p1, point p2) { } circle byThreePoints(point * warp) { - point center; - double radius, x, y; - double ang_a; - double ang_b; - - if (fabs(warp[1].x - warp[0].x) < exp || fabs(warp[2].x - warp[1].x) < exp) return (circle){.center=(point){0, 0}, .radius=-1}; - - ang_a = straightAngle(warp[1], warp[0]); - ang_b = straightAngle(warp[2], warp[1]); - - if (fabs(ang_b) < exp || fabs(ang_b - ang_a) < exp) { - return (circle){.center=(point){0, 0}, .radius=-1}; - } - - x = (ang_a * ang_b * (warp[0].y - warp[2].y) + ang_b * (warp[0].x + warp[1].x) - ang_a * (warp[1].x + warp[2].x)) / (2 * (ang_b - ang_a)); - y = (-(1/ang_b) * (x - (warp[1].x + warp[2].x) / 2) + (warp[1].y + warp[2].y) / 2); - center = (point){x, y}; + point center = getCenter(warp); + double radius; + center.x += warp[0].x, center.y += warp[0].y; radius = distance(center, warp[0]); + return (circle){center, radius}; } +point getCenter(point * warp) { + double sqrt1, sqrt2, scalar; + point vec1 = (point){warp[1].x - warp[0].x, warp[1].y - warp[0].y}; + point vec2 = (point){warp[2].x - warp[0].x, warp[2].y - warp[0].y}; + + sqrt1 = vec1.x * vec1.x + vec1.y * vec1.y; + sqrt2 = vec2.x * vec2.x + vec2.y * vec2.y; + scalar = (vec1.x * vec2.y - vec2.x * vec1.y) * 2; + + return (point){.x=((vec2.y * sqrt1 - vec1.y * sqrt2) / scalar), \ + .y=((vec1.x * sqrt2 - vec2.x * sqrt1) / scalar)}; +} + double straightAngle(point p1, point p2) { return (p1.y - p2.y) / (p1.x - p2.x); } @@ -85,7 +87,14 @@ double distance(point p1, point p2) { } void printCircle(circle crcl) { - printf("(x - %.4lf)^2 + (y - %.4lf)^2 = %.4lf^2\n", crcl.center.x, crcl.center.y, crcl.radius); - // printf("Center of circle at point (%.2lf, %.2lf)\nRadius is %.2lf\n", crcl.center.x, crcl.center.y, crcl.radius); + printf("(x - %.4lf)^2 + (y - %.4lf)^2 = %.4lf^2\n\n", crcl.center.x, crcl.center.y, crcl.radius); + printf("Center of circle at point (%.2lf, %.2lf)\nRadius is %.2lf\n", crcl.center.x, crcl.center.y, crcl.radius); } +int isCover(circle crcl, points pts) { + for (int i = 0; i < pts.length; ++i) { + if (!belongs(crcl, pts.array[i])) return i; + } + + return pts.length; +} \ No newline at end of file diff --git a/ComputationalGeometry/6Ex/SCP.h b/ComputationalGeometry/6Ex/SCP.h index 745251c..206675a 100644 --- a/ComputationalGeometry/6Ex/SCP.h +++ b/ComputationalGeometry/6Ex/SCP.h @@ -14,8 +14,10 @@ double powd(double number); circle primitive(point * N, int nlen); circle centermass(point p1, point p2); circle byThreePoints(point * warp); +point getCenter(point * warp); double straightAngle(point p1, point p2); double distance(point p1, point p2); void printCircle(circle crcl); +int isCover(circle crcl, points pts); #endif diff --git a/ComputationalGeometry/6Ex/hope.c b/ComputationalGeometry/6Ex/hope.c index 69eee1a..85b5c5c 100644 --- a/ComputationalGeometry/6Ex/hope.c +++ b/ComputationalGeometry/6Ex/hope.c @@ -2,49 +2,46 @@ circle hope(point * ps, int length) { circle crcl; + circle minimum = (circle){(point){0, 0}, 0}; + + if (length == 1) { + return (circle){ps[0], 0}; + } else if (length < 1) { + return minimum; + } for (int i = 0; i < length; ++i) { - point temp = ps[i]; - ps[i] = ps[length - 1]; + point temp = ps[0]; + + for (int i = 1; i < length; ++i) ps[i - 1] = ps[i]; ps[length - 1] = temp; for (int j = 0; j < length - 1; ++j) { - crcl = centermass(temp, ps[i]); - if (isSuit(crcl, ps, length - 1)) return crcl; + crcl = centermass(temp, ps[j]); + if (isCover(crcl, (points){ps, length}) == length) { + if (minimum.radius < exp || (minimum.radius - crcl.radius > exp)) { + minimum = crcl; + } + } } - - ps[length - 1] = ps[i]; - ps[i] = temp; } for (int i = 0; i < length; ++i) { - point temp = ps[i]; - ps[i] = ps[length - 1]; - ps[length - 1] = temp; - for (int j = 0; j < length - 1; ++j) { - temp = ps[j]; - ps[j] = ps[length - 2]; - ps[length - 2] = temp; - for (int k = 0; k < length - 2; ++k) { - point warp[] = {ps[length - 1], ps[length - 2], ps[k]}; + point warp[] = {ps[i], ps[j], ps[k]}; crcl = byThreePoints(warp); - if (fabs(crcl.radius + 1) > exp && isSuit(crcl, ps, length - 2)) return crcl; + if ((fabs(crcl.radius + 1) > exp) && (isCover(crcl, (points){ps, length}) == length)) { + if (minimum.radius < exp || (minimum.radius - crcl.radius > exp)) { + minimum = crcl; + } + } } - - temp = ps[j]; - ps[j] = ps[length - 2]; - ps[length - 2] = temp; } - - temp = ps[i]; - ps[i] = ps[length - 1]; - ps[length - 1] = temp; } - return (circle){.center=(point){0, 0}, .radius=0}; + return minimum; } bool isSuit(circle crcl, point * ps, int length) { diff --git a/ComputationalGeometry/6Ex/input.txt b/ComputationalGeometry/6Ex/input.txt index 871e879..8040f84 100644 --- a/ComputationalGeometry/6Ex/input.txt +++ b/ComputationalGeometry/6Ex/input.txt @@ -1 +1 @@ -1 0 -1 0 +36.8400 -33.0600 -11.2900 95.9600 -56.8500 -77.7900 61.9600 59.0600 -9.8400 -14.3700 1.0100 85.6400 \ No newline at end of file diff --git a/ComputationalGeometry/6Ex/main.c b/ComputationalGeometry/6Ex/main.c index b8c4fd0..fc5bdd7 100644 --- a/ComputationalGeometry/6Ex/main.c +++ b/ComputationalGeometry/6Ex/main.c @@ -6,6 +6,7 @@ int main(void) { points pts; point N[3]; circle crcl; + int numLose; pts = getPoints(); if (pts.array == NULL) return -1; @@ -16,9 +17,15 @@ int main(void) { crcl = MEC(pts.array, pts.length, N, 0); printCircle(crcl); + if ((numLose = isCover(crcl, pts)) == pts.length) printf("\nThe circle covers all the points!\n"); + else printf("\nThe circle misses a maximum of %d points", pts.length - numLose); + printf("\nReliable algorithm:\n"); crcl = hope(pts.array, pts.length); printCircle(crcl); + + if (isCover(crcl, pts) == pts.length) printf("\nThe circle covers all the points!\n"); + else printf("\nThe circle misses a maximum of %d points", pts.length - numLose); free(pts.array); diff --git a/Session/3Ex/ina.txt b/Session/3Ex/ina.txt new file mode 100644 index 0000000..9831405 --- /dev/null +++ b/Session/3Ex/ina.txt @@ -0,0 +1 @@ +1 3 7 8 11 \ No newline at end of file diff --git a/Session/3Ex/inb.txt b/Session/3Ex/inb.txt new file mode 100644 index 0000000..5292f80 --- /dev/null +++ b/Session/3Ex/inb.txt @@ -0,0 +1 @@ +-1 5 7 10 \ No newline at end of file diff --git a/Session/3Ex/main.c b/Session/3Ex/main.c new file mode 100644 index 0000000..23712c0 --- /dev/null +++ b/Session/3Ex/main.c @@ -0,0 +1,79 @@ +#include +#include +#include + +#define eps 1.e-6 + +double * getArray(FILE * file); +int level(double * arr1, int len1, double * arr2, int len2); +int contains(double * arr, int len, double num); + +double * getArray(FILE * file) { + int size = 2, i = 0; + double * arr = (double *)malloc(size * sizeof(double)); + double current; + + if (arr == NULL) return NULL; + + while (fscanf(file, "%lf", ¤t) == 1) { + if (size == (i + 1)) arr = (double *)realloc(arr, (size *= 2) * sizeof(double)); + arr[++i] = current; + } + + arr = (double *)realloc(arr, (i + 1) * sizeof(double)); + arr[0] = i; + + if (i == 0) return NULL; + return arr; +} + +int level(double * arr1, int len1, double * arr2, int len2) { + double start = 0, end = 0; + int count = 0, len = len1 < len2 ? len2 : len1; + + for (int i = 0; i < len; ++i) { + if (start - arr1[i % len2] > eps || (i == 0)) start = arr2[i % len2]; + if (arr1[i % len1] - end > eps || i == 0) end = arr1[i % len1]; + } + + + + for (int i = 0; i < len; ++i) { + if (i < len1) if (arr1[i] - start > -eps && end - arr1[i] > -eps && (contains(arr2, len2, arr1[i]) == 0)) ++count; + if (i < len2) if (arr2[i] - start > -eps && end - arr2[i] > -eps) ++count; + } + + return count; +} + +int contains(double * arr, int len, double num) { + for (int i = 0; i < len; ++i) if (fabs(arr[i] - num) < eps) return 1; + return 0; +} + +int main(void) { + FILE * file = fopen("ina.txt", "r"); + double *arr1, *arr2; + + if (file == NULL) return -1; + arr1 = getArray(file); + if (fclose(file) != 0) return -1; + if (arr1 == NULL) return -1; + + file = fopen("inb.txt", "r"); + if (file == NULL) return -1; + arr2 = getArray(file); + if (fclose(file) != 0) return -1; + if (arr1 == NULL) return -1; + + file = fopen("output.txt", "w"); + if (file == NULL) return -1; + if (fprintf(file, "%d", level(&arr1[1], (int)arr1[0], &arr2[1], (int)arr2[0])) <= 0) {fclose(file); return -1;} + if (fclose(file) != 0) return -1; + + free(arr1); + free(arr2); + + return 0; +} + diff --git a/Session/3Ex/makefile b/Session/3Ex/makefile new file mode 100644 index 0000000..f0ffbed --- /dev/null +++ b/Session/3Ex/makefile @@ -0,0 +1,37 @@ +CFLAGS = -mfpmath=sse \ + -fstack-protector-all \ + -W \ + -Wall \ + -Wextra \ + -Wunused \ + -Wcast-align \ + -Werror \ + -pedantic \ + -pedantic-errors \ + -Wfloat-equal \ + -Wpointer-arith \ + -Wformat-security \ + -Wmissing-format-attribute \ + -Wformat=1 \ + -Wwrite-strings \ + -Wcast-align \ + -Wno-long-long \ + -std=gnu99 \ + -Wstrict-prototypes \ + -Wmissing-prototypes \ + -Wmissing-declarations \ + -Wold-style-definition \ + -Wdeclaration-after-statement \ + -Wbad-function-cast \ + -Wnested-externs \ + -O3 \ + -D_DEBUG -g \ + -c + +all: main.o + gcc main.o -lssp && del *.o + a.exe + +main.o: main.c + gcc $(CFLAGS) main.c + diff --git a/Session/4Ex/ina.txt b/Session/4Ex/ina.txt new file mode 100644 index 0000000..27290a3 --- /dev/null +++ b/Session/4Ex/ina.txt @@ -0,0 +1 @@ +89 12 65 \ No newline at end of file diff --git a/Session/4Ex/inb.txt b/Session/4Ex/inb.txt new file mode 100644 index 0000000..517e4c0 --- /dev/null +++ b/Session/4Ex/inb.txt @@ -0,0 +1 @@ +5 4 3 2 1 0 \ No newline at end of file diff --git a/Session/4Ex/main.c b/Session/4Ex/main.c new file mode 100644 index 0000000..34fb16c --- /dev/null +++ b/Session/4Ex/main.c @@ -0,0 +1,93 @@ +#include +#include +#include + +#define eps 1.e-6 + +int ctv(const void * unum1, const void * unum2); +double * getArray(FILE * file); +int isContains(double * arr1, int len1, double * arr2, int len2); +int write(const char * str); + +int ctv(const void * unum1, const void * unum2) { + double num1 = *(double *)unum1; + double num2 = *(double *)unum2; + + if (num1 - num2 > eps) return 1; + else if (num2 - num1 > eps) return -1; + else return 0; +} + +double * getArray(FILE * file) { + int size = 2, i = 0; + double * arr = (double *)malloc(size * sizeof(double)); + double current; + + if (arr == NULL) return NULL; + + + while (fscanf(file, "%lf", ¤t) == 1) { + if (size == (i + 1)) { + arr = (double *)realloc(arr, (size *= 2) * sizeof(double)); + if (arr == NULL) return NULL; + } + arr[++i] = current; + } + + arr[0] = i; + return arr; +} + +int isContains(double * arr1, int len1, double * arr2, int len2) { + for (int i = 0, k = 0; i < len2; ++i, ++k) { + if (fabs(arr2[i] - arr1[k]) < eps) { + if (k + 1 == len1) return 1; + } else if (len2 - i <= len1) return 0; + else k = -1; + } + + return 0; +} + +int write(const char * str) { + FILE * file = fopen("output.txt", "w"); + + if (file == NULL) return -1; + if (fprintf(file, "%s", str) < 1) {fclose(file); return -1;} + if (fclose(file) != 0) return -1; + + return 0; +} + +int main(void) { + FILE * file = fopen("ina.txt", "r"); + double *arr1, *arr2; + + if (file == NULL) return -1; + arr1 = getArray(file); + if (arr1 == NULL) {fclose(file); return -1;} + if (fclose(file) != 0) return -1; + + file = fopen("inb.txt", "r"); + if (file == NULL) return -1; + arr2 = getArray(file); + if (arr2 == NULL) {fclose(file); return -1;} + if (fclose(file) != 0) return -1; + + if ((int)arr2[0] < (int)arr1[0]) return -1; + + qsort(&arr1[1], (int)arr1[0], sizeof(double), ctv); + qsort(&arr2[1], (int)arr2[0], sizeof(double), ctv); + + if (isContains(&arr1[1], (int)arr1[0], &arr2[1], (int)arr2[0]) == 0) { + free(arr1); + free(arr2); + + return write("NO"); + } else { + free(arr1); + free(arr2); + + return write("YES"); + } +} diff --git a/Session/4Ex/makefile b/Session/4Ex/makefile new file mode 100644 index 0000000..f0ffbed --- /dev/null +++ b/Session/4Ex/makefile @@ -0,0 +1,37 @@ +CFLAGS = -mfpmath=sse \ + -fstack-protector-all \ + -W \ + -Wall \ + -Wextra \ + -Wunused \ + -Wcast-align \ + -Werror \ + -pedantic \ + -pedantic-errors \ + -Wfloat-equal \ + -Wpointer-arith \ + -Wformat-security \ + -Wmissing-format-attribute \ + -Wformat=1 \ + -Wwrite-strings \ + -Wcast-align \ + -Wno-long-long \ + -std=gnu99 \ + -Wstrict-prototypes \ + -Wmissing-prototypes \ + -Wmissing-declarations \ + -Wold-style-definition \ + -Wdeclaration-after-statement \ + -Wbad-function-cast \ + -Wnested-externs \ + -O3 \ + -D_DEBUG -g \ + -c + +all: main.o + gcc main.o -lssp && del *.o + a.exe + +main.o: main.c + gcc $(CFLAGS) main.c + diff --git a/Session/5Ex/main.c b/Session/5Ex/main.c new file mode 100644 index 0000000..8139c7f --- /dev/null +++ b/Session/5Ex/main.c @@ -0,0 +1,70 @@ +#include +#include + +int cmp(const * void unum1, const void * unum2) { + double num1 = *(double *)unum1; + double num2 = *(double *)unum2; + + if (num1 - num2 > eps) return -1; + else if (num2 - num1 > eps) return 1; + + return 0; +} + +double * getArray(FILE * file) { + int size = 2, i = 0; + double * arr = (double *)malloc(size * sizeof(double)); + double current; + + if (arr == NULL) return NULL; + + while (fscanf(file, "%lf", ¤t) == 1) { + if (size == i + 1) { + arr = (double *)realloc(arr, (size *= 2) * sizeof(double)); + if (arr == NULL) return NULL; + } + + arr[++i] = current; + } + + arr[0] = i; + return arr; +} + +int * compare(double *arr1, int len1, double *arr2, int len2) { + int *answer = (int *)malloc(2 * sizeof(int)); + int len = len1 < len2 ? len1 : len2; + + if (answer == NULL) return NULL; + answer[0] = 0, answer[1] = 0; + + for (int i = 0; i < len; ++i) { + + } + + return answer; +} + +int main(void) { + FILE * file = fopen("ina.txt", "r"); + double *arr1, *arr2; + + if (file == NULL) return -1; + arr1 = getArray(file); + if (arr1 == NULL) {fclose(file); return -1;} + if (fclose(file) != 0) return -1; + + file = fopen("inb.txt", "r"); + if (file == NULL) return -1; + arr2 = getArray(file); + if (arr2 == NULL) {fclose(file); return -1;} + if (Fclose(file) != 0) return -1; + + qsort(&arr1[1], (int)arr1[0], sizeof(double), cmp); + qsort(&arr2[1], (int)arr2[0], sizeof(double), cmp); + + free(arr1); + free(arr2); + + return 0; +}