Численное интегрирование

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
 
#ifndef NDEBUG
#define DEBUG(EXPR, ...)                        \
    fprintf(stderr, EXPR "\n", __VA_ARGS__)     
#else
#define DEBUG(EXPR, ...) ((void) 0)
#endif
 
double square(double x)
{
    return x * x;
}
 
#define FUNCTION_TABLE                          \
    DEF_FUNC(sin, "sin x")                      \
    DEF_FUNC(cos, "cos x")                      \
    DEF_FUNC(tan, "tan x")                      \
    DEF_FUNC(exp, "e^x")                        \
    DEF_FUNC(square, "x^2")                     
 
typedef double (*FUNC)(double);
 
#define DEF_FUNC(fun, name) fun,
static FUNC functions[] = 
{
    FUNCTION_TABLE
};
#undef DEF_FUNC
 
#define DEF_FUNC(fun, name) name,
static const char *names[] =
{
    FUNCTION_TABLE
};
#undef DEF_FUNC
 
#define EPS 5e-4
#define FUNCTIONS_COUNT (sizeof functions / sizeof *functions)
#define INITIAL_STEP 1e-2
#define STEP_TOLERANCE 10
#define THETA 0.333             /* For trapezoid */
 
double trapezoid(FUNC func, double x, double step)
{
    return (func(x) + func(x + step)) * step / 2.0;
}
 
double runge_estimation(double i2, double i)
{
    return THETA * fabs(i2 - i);
}
 
enum action
{
    INCREASE, DECREASE
};
 
void adjust_step(double *current_step, enum action action)
{
    *current_step *= action == DECREASE ? .5 : 2;
}
 
double integrate(FUNC func, double a, double b)
{
    double result = 0;
    double step = INITIAL_STEP, step_2;
    double x = a;
    
    size_t num_steps = 0;
    
    while(x < b)
    {
        if(x + step > b)
            step = b - x;
 
        step_2 = step / 2;
        
        double integr = trapezoid(func, x, step);
        double integr_2 = trapezoid(func, x, step_2) +
            trapezoid(func, x + step_2, step_2);
 
        double error_estimation = runge_estimation(integr, integr_2);
        double acceptable_error = step * EPS / (b - a);
        
        if(acceptable_error < error_estimation) /* Step is too large */
            adjust_step(&step, DECREASE);
        else
        {
            result += integr_2;
            x += step;
            ++num_steps;
 
            /* Step is too small */
            if(acceptable_error >= STEP_TOLERANCE * error_estimation)
                adjust_step(&step, INCREASE);
        }
    }
 
    DEBUG("Integral = %.3f (%zu steps)", result, num_steps);
    
    return result;
}
 
int menu(void)
{
    size_t i;
    size_t num;
    
    puts("Select a function to integrate");
    
    for(i = 0; i < FUNCTIONS_COUNT; ++i)
        printf("% 4d> %s\n", (int) i + 1, names[i]);
    puts("or enter 0 to exit");
 
    printf("\n>>> ");
    fflush(stdout);
 
    if(scanf("%zu", &num) != 1)
        return -1;
 
    return num;
}
 
int main(void)
{
    int choice = menu();
 
    if(choice == -1)
    {
        fputs("Incorrect input\n", stderr);
        exit(1);
    }
 
    if(choice > (int) FUNCTIONS_COUNT)
    {
        fprintf(stderr, "Your choice must be in the range [0, %zu]\n", FUNCTIONS_COUNT);
        exit(1);
    }
 
    if(choice != 0)
    {
        double a, b;
 
        fputs("Left bound: ", stdout);
        fflush(stdout);
        scanf("%lf", &a);
 
        fputs("Right bound: ", stdout);
        fflush(stdout);
        scanf("%lf", &b);
                
        double result = integrate(functions[choice - 1], a, b);
 
        printf("Integral of f(x) = %s in range [%.2f, %.2f] is %.3f\n",
               names[choice - 1], a, b, result);
    }
    
    exit(0);
}

#define false 0
#define true 1
typedef int bool;