Math.lab3

#include <stdio.h>
#include <math.h>
#define F(x) (x)*(x)*(x)*(x)*(x) - 3*(x)*(x) + 1
#define f(x) 5*(x)*(x)*(x)*(x) - 6*(x)
#define f_2(x) 20*(x)*(x)*(x) - 6
#define EPS_BIN 0.1
#define EPS 0.0005

float half_division(float, float);
float newton(float);
void kantorovich(float);

int main()
{
    float a0 = half_division(-1.0, 0.0);
    float a1 = half_division(0.0, 1.0);
    float a2 = half_division(1.0, 2.0);

    printf("\nRoots approximations (eps = 0.1):\n");
    printf("\tx1 = %f\tF(x1) = %f\n", a0, F(a0));
    kantorovich(a0);
    printf("\tx2 = %f\tF(x2) = %f\n", a1, F(a1));
    kantorovich(a1);
    printf("\tx3 = %f\tF(x3) = %f\n", a2, F(a2));
    kantorovich(a2);
    putchar('\n');

    float x1 = newton(a0);
    float x2 = newton(a1);
    float x3 = newton(a2);
    printf("Roots approximations (eps = 0.0005):\n");
    printf("\tx1 = %f\tF(x1) = %f\n\n", x1, F(x1));
    printf("\tx2 = %f\tF(x2) = %f\n\n", x2, F(x2));
    printf("\tx3 = %f\tF(x3) = %f\n\n", x3, F(x3));

    return 0;
}

void kantorovich(float x0)
{
    float fx0 = f(x0);
    float B = fabs(1/fx0);
    float Fx0 = F(x0);
    float nu = fabs(Fx0/fx0);
    float f_2 = f_2(x0);
    float K = fabs(f_2);
    float h = B * nu * K;
    float a = (1 - sqrt(1 - 2*h)) / h * nu;
    if ((h < 0.5) && (a < 0.5))
        printf("\th = %f < 0.5\ta(h) = %f < 0.5\n\n\n", h, a);
}

float newton(float x0)
{
    float fx0 = f(x0);
    float Fx0 = F(x0);
    float x1 = x0 - Fx0 / fx0;

    while (fabs(x0 - x1) >= EPS) {
        x0 = x1;
        Fx0 = F(x0);
        x1 = x0 - Fx0 / fx0;
    }

    return x1;
}

float half_division(float a, float b)
{
    float left = a, right = b;
    float middle = left + (right-left)/2.0;

    while ((right-left) > EPS_BIN/2.0){
        if (F(middle) < 0) {
            if (F(left) > 0)
                right = middle;
            else
                left = middle;
        } else {
            if (F(left) < 0)
                right = middle;
            else
                left = middle;
        }
        middle = left + (right-left)/2.0;
    }

    return middle;
}