diff --git a/2025.05.02/06Ex/Makefile b/2025.05.02/06Ex/Makefile index c209389..b4b5f65 100644 --- a/2025.05.02/06Ex/Makefile +++ b/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 $@ diff --git a/2025.05.02/06Ex/array_io.c b/2025.05.02/06Ex/array_io.c new file mode 100644 index 0000000..0ca64d1 --- /dev/null +++ b/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; +} diff --git a/2025.05.02/06Ex/array_io.h b/2025.05.02/06Ex/array_io.h new file mode 100644 index 0000000..15cebe3 --- /dev/null +++ b/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 diff --git a/2025.05.02/06Ex/comp.h b/2025.05.02/06Ex/comp.h new file mode 100644 index 0000000..6b09d43 --- /dev/null +++ b/2025.05.02/06Ex/comp.h @@ -0,0 +1,50 @@ +#ifndef COMP_H +#define COMP_H + +#include +#include + +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 diff --git a/2025.05.02/06Ex/main.c b/2025.05.02/06Ex/main.c index 442f122..11c7d05 100644 --- a/2025.05.02/06Ex/main.c +++ b/2025.05.02/06Ex/main.c @@ -1,48 +1,65 @@ #include +#include #include #include #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; - int m, k, cl, it, task = 1; + double t = 0, *d, a, b, eps, x = 0; + 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) && - sscanf(argv[1], "%lf", &a) == 1 && - sscanf(argv[2], "%lf", &b) == 1 && + !((argc == 7) && + sscanf(argv[1], "%d", °_poly) == 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], "%d", &m) == 1) && m > 0) && - ((sscanf(argv[5], "%d", &k) == 1) && ((0 <= k) && (k <= len_f)))) + (sscanf(argv[4], "%lf", &eps) == 1 && (eps >= 0)) && + ((sscanf(argv[5], "%d", &num_iter) == 1) && num_iter > 0) && + ((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]; + + it = init_array(f, d, deg_poly, a, b, eps, &x); + + if (!it) + { + t = clock(); + it = t6_solve(f, deg_poly, d, a, b, eps, num_iter, &x); + t = (clock() - t) / CLOCKS_PER_SEC; + } - t = clock(); - it = t1_solve(f, a, b, eps, m, &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; } } diff --git a/2025.05.02/06Ex/polynom.c b/2025.05.02/06Ex/polynom.c new file mode 100644 index 0000000..3a68254 --- /dev/null +++ b/2025.05.02/06Ex/polynom.c @@ -0,0 +1,43 @@ +#include "polynom.h" + +#include +#include +#include + +// 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; +} + diff --git a/2025.05.02/06Ex/polynom.h b/2025.05.02/06Ex/polynom.h new file mode 100644 index 0000000..65611c0 --- /dev/null +++ b/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 diff --git a/2025.05.02/06Ex/solve.c b/2025.05.02/06Ex/solve.c index 233ee33..b16f469 100644 --- a/2025.05.02/06Ex/solve.c +++ b/2025.05.02/06Ex/solve.c @@ -1,103 +1,57 @@ #include "solve.h" +#include "comp.h" +#include "polynom.h" #include #include -#include -#include -#include - -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 ) { - const int len = m + 1; + 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; - double last_y = f(a); - 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) + int it; + for (it = 1; it <= M; ++it) { - 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; - x_lst[i] = t_x; + if (is_eps(y, eps)) + { + *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; } + diff --git a/2025.05.02/06Ex/solve.h b/2025.05.02/06Ex/solve.h index 9b10deb..5042f7a 100644 --- a/2025.05.02/06Ex/solve.h +++ b/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