david

#include <stdio.h>
#include <math.h>

double fST(double s, double t){
    return 1/sqrt(pow(s, 4) + pow(t, 4));
}

double integrate_dt(double s, double a, double b, double n){
    double h = (b-a)/n;
    double sum = 0;

    for (int i = 1; i <= n / 2; i++)
        sum += 4 * fST(s, a + (2 * i - 1) * h);

    for (int i = 1; i < n / 2; i++)
        sum += 2 * fST(s, a + 2 * i * h);

    return (fST(s, a) + sum + fST(s, b)) * h / 3;
}
double get_fS(double s, double a, double b, long n, double eps) {
    double x = integrate_dt(s, a, b, n);
    double y = integrate_dt(s, a, b, 2 * n);

    if (fabs(x - y) < eps)
        return y;

    return get_fS(s, a, (a+b)/2, n, eps) + get_fS(s, (a + b) / 2, b, n, eps);
}
double fS(double s, double eps){
    return get_fS(s, 0, 1, 2, eps);
}
//check fS value when s is some const here
//сheking succseed

//now we have some f(s) func
//the problem is now to find x for integral f(s) ds from 1 to x = alpha

double integrate(double a, double b, long n, double eps){
    double h = (b-a)/n;
    double sum = 0;

    for (int i = 1; i <= n / 2; i++)
        sum += 4 * fS(a + (2 * i - 1) * h, eps);

    for (int i = 1; i < n / 2; i++)
        sum += 2 * fS(a + 2 * i * h, eps);

    return (fS(a, eps) + sum + fS(b, eps)) * h / 3;
}
double getIntegral(double a, double b, long n, double eps) {
    double x = integrate(a, b, n, eps);
    double y = integrate(a, b, 2 * n, eps);

    if (fabs(x - y) < eps)
        return y;
    //printf("!\n");
    return getIntegral(a, (a+b)/2, n, eps) + getIntegral((a + b) / 2, b, n, eps);
}
double f(double x, double alpha, double eps){
    return (getIntegral(1, x, 2, eps) - alpha);
}
//check f + alpha value
//checking succeed

//f is now a func we need an x for
//so lets make a dev for it and then use the newton method

double df(double x, double alpha, double eps){
    return (f(x+eps, alpha, eps) - f(x-eps, alpha, eps))/(2*eps);
}
//check df value

double newtone_method(double x0, double alpha, double eps) {
    double x;
    do {
        x = x0;
        x0 -= f(x0, alpha, eps) / df(x0, alpha, eps);
    } while (fabs(x - x0) > eps);
    return x0;
}

int main(){
    double alpha;
    double eps;
    printf("input alpha: ");
    scanf("%lf", &alpha);
    if(alpha >=1 || alpha <0){
        printf("invalid alpha: integral cant be 1 or more or less than 0\n");
        return -1;
    }
    printf("input eps: ");
    scanf("%lf", &eps);
    if(eps > 1 || eps < 1e-8){
        printf("invalid eps\n");
        return -1;
    }
    double x = newtone_method(1, alpha, eps);
    printf("x = %lf\nf(x) = %lf", x, f(x+eps, alpha, eps));

    //checking fS
    //printf("%lf", fS(1, eps));
    //UPD succeed

    //checking f
    //printf("%lf\n", f(1+eps, alpha, eps));
    //UPD succeed
    //cheking df
    //printf("%lf\n", df(2, alpha, eps));
    //printf("%lf\n", fS(2, eps));
    //UPD succeed
    return 0;
}