Написал надёжный алгоритм поиска наименьшего круга

ComputationalGeometry/6Ex/SCP.c

@@ -56,8 +56,17 @@ circle centermass(point p1, point p2) {
 circle byThreePoints(point * warp) {
 	point center;
 	double radius, x, y;
-	double ang_a = straightAngle(warp[1], warp[0]);
-	double ang_b = straightAngle(warp[2], warp[1]);
+	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);

ComputationalGeometry/6Ex/hope.c

@@ -0,0 +1,53 @@
+#include "hope.h"
+
+circle hope(point * ps, int length) {
+	circle crcl;
+	
+	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) {
+			crcl = centermass(temp, ps[i]);
+			if (isSuit(crcl, ps, length - 1)) return 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]};
+				crcl = byThreePoints(warp);
+
+				if (fabs(crcl.radius + 1) > exp && isSuit(crcl, ps, length - 2)) return 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};
+}
+
+bool isSuit(circle crcl, point * ps, int length) {
+	for (int i = 0; i < length; ++i) if (!belongs(crcl, ps[i])) return false;
+	return true;
+}

ComputationalGeometry/6Ex/hope.h

@@ -0,0 +1,10 @@
+#ifndef HOPE
+#define HOPE
+
+#include <stdbool.h>
+#include "SCP.h"
+
+circle hope(point * ps, int length);
+bool isSuit(circle crcl, point * ps, int length);
+
+#endif

ComputationalGeometry/6Ex/main.c

@@ -1,5 +1,6 @@
 #include "tools.h"
 #include "SCP.h"
+#include "hope.h"
 
 int main(void) {
 	points ps;
@@ -11,11 +12,19 @@ int main(void) {
 	ps = getPoints(file);
 	if (ps.array == NULL) return -2;
 	
+	printf("\nPoints:\n");	
 	for (int i = 0; i < ps.length; ++i) printf("(%.2lf, %.2lf) ", ps.array[i].x, ps.array[i].y);
-	printf("\n");
-
+	printf("\n\n");
+	
+	printf("Fast algorithm:\n");
 	crcl = MEC(ps.array, ps.length, N, 0);
 	printCircle(crcl);
+	
+	printf("\n");
+	
+	printf("Reliable algorithm:\n");
+	crcl = hope(ps.array, ps.length);
+	printCircle(crcl);
 
 	free(ps.array);	
 

ComputationalGeometry/6Ex/makefile

@@ -28,8 +28,8 @@ CFLAGS = -mfpmath=sse \
 		 -D_DEBUG -g \
 		 -c
 
-all: main.o SCP.o tools.o
-	gcc main.o SCP.o tools.o -lssp && del *.o
+all: main.o SCP.o hope.o tools.o
+	gcc main.o SCP.o hope.o tools.o -lssp && del *.o
 	a.exe
 
 main.o: main.c
@@ -38,5 +38,8 @@ main.o: main.c
 SCP.o: SCP.c
 	gcc $(CFLAGS) SCP.c
 
+hope.o: hope.c
+	gcc $(CFLAGS) hope.c
+
 tools.o: tools.c
 	gcc $(CFLAGS) tools.c