Started task 6

2025.05.02/06Ex/Makefile

@@ -21,7 +21,7 @@ else
 endif
 
 TARGET = a06.$(EXE)
-OBJ = main.o solve.o init_f.o
+OBJ = main.o solve.o init_f.o array_io.o polynom.o
 
 %.o: %.c
 	gcc $(WFLAGS) $(LDFLAGS) -c $< -o $@

2025.05.02/06Ex/array_io.c

@@ -0,0 +1,41 @@
+#include "array_io.h"
+#include "comp.h"
+
+int init_array (
+		double (*f) (double), 
+		double *d, int m,
+		double a, double b,
+	       	double eps, double *ret
+) {
+	const double h = (b - a) / m;
+	const int len = m+1;
+
+	double *y_lst = d;
+	double *x_lst = d + len;
+	double *t_lst = d + (len < 1);
+
+	double x = a;
+
+	for (int i = 0; i < len; ++i)
+	{
+		double y = f(x);
+
+		if (is_eps(y, eps))
+		{
+			*ret = x;
+			return 1;
+		}
+
+		for (int j = 0; j < i; ++j)
+			if (is_equal(y, y_lst[j]))
+				return -1;
+
+		y_lst[i] = y;
+		x_lst[i] = x;
+		t_lst[i] = x;
+
+		x += h;
+	}
+
+	return 0;
+}

2025.05.02/06Ex/array_io.h

@@ -0,0 +1,11 @@
+#ifndef ARRAY_IO_H
+#define ARRAY_IO_H
+
+int init_array (
+		double (*f) (double), 
+		double *d, int m,
+		double a, double b,
+	       	double eps, double *ret
+);
+
+#endif

2025.05.02/06Ex/comp.h

@@ -0,0 +1,50 @@
+#ifndef COMP_H
+#define COMP_H
+
+#include <float.h>
+#include <math.h>
+
+static inline double * fpmax (double *pa, double *pb, double fa, double fb, double *max_f_p)
+{
+	if ((fa - fb) > DBL_EPSILON) 
+	{
+		*max_f_p = fa;
+		return pa;
+	} else
+	{
+		*max_f_p = fb;
+		return pb;
+	}
+}
+
+static inline double * fp_abs_max (double *pa, double *pb, double *fa, double *fb, double **max_f_p)
+{
+	if ((fabs(*fa) - fabs(*fb)) > DBL_EPSILON) 
+	{
+		*max_f_p = fa;
+		return pa;
+	} else
+	{
+		*max_f_p = fb;
+		return pb;
+	}
+}
+
+static inline int is_equal (const double a, const double b)
+{
+	double diff = a - b;
+	double max_val = (a > b) ? a : b;
+	return ((diff < 0) ? -diff : diff) < (DBL_EPSILON * max_val);
+}
+
+static inline int is_null (const double a)
+{
+	return ((a < 0) ? -a : a) < DBL_EPSILON;
+}
+
+static inline int is_eps (const double a, const double eps)
+{
+	return (((a < 0) ? -a : a) - eps) < DBL_EPSILON;
+}
+
+#endif

2025.05.02/06Ex/main.c

@@ -1,48 +1,65 @@
 #include <stdio.h>
+#include <stdlib.h>
 #include <time.h>
 #include <math.h>
 
 #include "init_f.h"
+#include "array_io.h"
 #include "solve.h"
 
-/* ./a.out a b eps M k */
+/* ./a.out m a b eps M k */
 int main(int argc, char *argv[])
 {
-	double t, a, b, eps, x = 0;
+	double t = 0, *d, a, b, eps, x = 0;
-	int m, k, cl, it, task = 1;
+	int deg_poly, num_iter, k, cl, it, task = 6;
 	
 	double (*f) (double);
 	double (*f_lst[]) (double) = {f0, f1, f2, f3, f4, f5, f6};
 	int len_f = sizeof(f_lst) / sizeof(f_lst[0]);
 
 	if (
-			!((argc == 6) && 
+			!((argc == 7) && 
-			sscanf(argv[1], "%lf", &a) == 1 &&
+			sscanf(argv[1], "%d", &deg_poly) == 1 &&
-			sscanf(argv[2], "%lf", &b) == 1 &&
+			sscanf(argv[2], "%lf", &a) == 1 &&
+			sscanf(argv[3], "%lf", &b) == 1 &&
 			(a <= b) &&
-			(sscanf(argv[3], "%lf", &eps) == 1 && (eps >= 0)) &&
+			(sscanf(argv[4], "%lf", &eps) == 1 && (eps >= 0)) &&
-			((sscanf(argv[4], "%d", &m) == 1) && m > 0) &&
+			((sscanf(argv[5], "%d", &num_iter) == 1) && num_iter > 0) &&
-			((sscanf(argv[5], "%d", &k) == 1) && ((0 <= k) && (k <= len_f))))
+			((sscanf(argv[6], "%d", &k) == 1) && ((0 <= k) && (k <= len_f))))
 	) {
-		fprintf(stderr, "Usage: %s a b eps M k\n", argv[0]);
+		fprintf(stderr, "Usage: %s m a b eps M k\n", argv[0]);
 		return -1;
 	}
 
+	d = (double *)malloc(3 * (deg_poly + 1) * sizeof(double));
+	if (!d)
+	{
+		fprintf(stderr, "Error: Not enough memory!\n");
+		return -2;
+	}
+	
 	f = f_lst[k];
 	
-	t = clock();
+	it = init_array(f, d, deg_poly, a, b, eps, &x);
-	it = t1_solve(f, a, b, eps, m, &x);
+	
-	t = (clock() - t) / CLOCKS_PER_SEC;
+	if (!it)
+	{
+		t = clock();
+		it = t6_solve(f, deg_poly, d, a, b, eps, num_iter, &x);
+		t = (clock() - t) / CLOCKS_PER_SEC;
+	}
 
 	cl = get_call_count();
 	
+	free(d);
+
 	if (it < 0)
 	{
 		fprintf(stdout, "%s : Task = %d NOT FOUND Count = %d T = %.2f\n", argv[0], task, cl, t);
 		return -2;
 	} else 
 	{
-		printf("%s : Task = %d X = %e Res = %e Its = %d Count = %d T = %.2f\n", argv[0], task, x, f(x), it, cl, t);
+		fprintf(stdout, "%s : Task = %d X = %e Res = %e Its = %d Count = %d T = %.2f\n", argv[0], task, x, f(x), it, cl, t);
 		return 0;
 	}
 }

2025.05.02/06Ex/polynom.c

@@ -0,0 +1,43 @@
+#include "polynom.h"
+
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+
+// the Newton interpolation polynomial
+double construct_poly (const double x_0, const int n, const double * restrict X, double * restrict Y)
+{
+	double value, start_value;
+
+	for (int k = 0; k < n-1; ++k)
+	{
+		double last_x; 
+		double last_y = Y[n-1];
+		
+		for (int i = n-2; i >= k; --i)
+		{
+			const double x_i = X[i-k];
+			const double y_i = Y[i];
+			last_x = X[i+1];
+			
+			if (fabs(last_x - x_i) < DBL_EPSILON)
+				return DBL_MAX;
+			
+			Y[i+1] = (last_y - y_i) / (last_x - x_i);
+
+			last_y = y_i;
+		}
+	}
+
+	start_value = 1;
+	value = 0;
+
+	for (int i = 0; i < n; ++i)
+	{
+		value += Y[i] * start_value;
+		start_value *= (x_0 - X[i]);
+	}
+
+	return value;
+}
+

2025.05.02/06Ex/polynom.h

@@ -0,0 +1,6 @@
+#ifndef POLYNOM_H
+#define POLYNOM_H
+
+double construct_poly (const double x_0, const int n, const double * restrict X, double * restrict Y);
+
+#endif

2025.05.02/06Ex/solve.c

@@ -1,103 +1,57 @@
 #include "solve.h"
+#include "comp.h"
+#include "polynom.h"
 
 #include <math.h>
 #include <float.h>
-#include <stdbool.h>
-#include <stdint.h>
-#include <string.h>
-
-int t1_solve (
-		double (*f) (double), 
-		double a, double b, 
-		double eps, int m, double *x
-) {
-	int it = 0;
-	uint64_t bits;
-	double c = DBL_MAX, y, y_a = f(a), y_b = f(b);
-	bool sgn_a, sgn_b, sgn_c;
-
-	memcpy(&bits, &y_a, sizeof(bits));
-	sgn_a = (bits >> 63) & 1;
-	memcpy(&bits, &y_b, sizeof(bits));
-	sgn_b = (bits >> 63) & 1;
-
-	if (fabs(y_a) - eps < DBL_EPSILON)
-	{
-		*x = a;
-		return 1;
-	} if (fabs(y_b) - eps < DBL_EPSILON)
-	{
-		*x = b;
-		return 1;
-	}
-
-	if (sgn_a == sgn_b)
-	{
-		*x = DBL_MAX;
-		return -1;
-	}
-	
-	for (it = 1; it <= m; ++it)
-	{
-		c = (a + b) * 0.5;
-		y = f(c);
-		
-		memcpy(&bits, &y, sizeof(bits));
-		sgn_c = (bits >> 63) & 1;
-		
-		if (fabs(y) - eps < DBL_EPSILON)
-			break;
-		else if ((fabs(c - a) < DBL_EPSILON) || (fabs(c - b) < DBL_EPSILON))
-			it = m+1;
-		else if (sgn_c == sgn_a)
-		{
-			a = c;
-			y_a = y;
-		} else if (sgn_c == sgn_b)
-		{
-			b = c;
-			y_b = y;
-		}
-	}
-	
-	if (it > m)
-		it = -1;
-	
-	*x = c;
-
-	return it;
-}
-
-
 
 int t6_solve (
 		double (*f) (double), 
-		const int m, double *d
+		const int m, double *d,
 		double a, double b, 
-		const double eps, const int M, double *x
+		const double eps, const int M, double *res
 ) {
+	double *maximum = d;
+
 	const int len = m + 1;
 	
 	double *y_lst = d;
 	double *x_lst = d + len;
 	double *t_lst = d + (len << 1);
 
-	double last_x = a;
+	int it;
-	double last_y = f(a);
+	for (it = 1; it <= M; ++it)
-	const int h = (b - a) / m;
-	
-	y_lst[0] = last_y;
-	x_lst[0] = last_x;
-	t_lst[0] = last_x;
-
-	for  (int i = 1, double t_x = a ; i < len ; i++, t_x += h)
 	{
-		double t_y = f(t_x);
+		double x = construct_poly(0, len, y_lst, x_lst);
+		double y = f(x);
 
-		y_lst[i] = t_y;
+		if (is_eps(y, eps))
-		x_lst[i] = t_x;
+		{
+			*res = x;
+			return it;
+		}
 
+		// Можно возвращение значений функции в y_lst можно было встроить в суммирование полинома, но мало толку
+		for (int i = 0; i < len; ++i)
+		{
+			double yi = t_lst[i];
+			y_lst[i] = yi;
 
+			if (is_equal(y, yi))
+				return -1;
+			
+			if ((fabs(yi) - fabs(*maximum)) > DBL_EPSILON)
+				maximum = &yi;
+		}
+
+		*maximum = y;
+		*(maximum + len) = x;
+		*(maximum + (len << 1)) = y;
 	}
 
+	(void)a;
+	(void)b;
+
+	return -1;
 }
+

2025.05.02/06Ex/solve.h

@@ -1,10 +1,11 @@
 #ifndef SOLVE_H
 #define SOLVE_H
 
-int t1_solve (
+int t6_solve (
 		double (*f) (double), 
+		const int m, double *d,
 		double a, double b, 
-		double eps, int m, double *x
+		const double eps, const int M, double *res
-		);
+);
 
 #endif