double integral(a, b, f, eps)
double a, b; /* концы отрезка */
double eps; /* требуемая точность */
double (*f)(); /* подынтегральная функция */
{
register long i;
double fab = (*f)(a) + (*f)(; /* сумма на краях */
double h, h2; /* шаг и удвоенный шаг */
long n, n2; /* число точек разбиения и оно же удвоенное */
double Sodd, Seven; /* сумма значений f в нечетных и в
четных точках */
double S, Sprev;/* значение интеграла на данной
и на предыдущей итерациях */
double x; /* текущая абсцисса */

niter = 0;
n = 10L; /* четное число */
n2 = n * 2;

h = fabs(b - a) / n2; h2 = h * 2.0;

/* Вычисляем первое приближение */
/* Сумма по нечетным точкам: */
for( Sodd = 0.0, x = a+h, i = 0;
i < n;
i++, x += h2 )
Sodd += (*f)(x);

/* Сумма по четным точкам: */
for( Seven = 0.0, x = a+h2, i = 0;
i < n-1;
i++, x += h2 )
Seven += f(x);

/* Предварительное значение интеграла: */
S = h / 3.0 * (fab + 4.0 * Sodd + 2.0 * Seven );
do{
niter++;
Sprev = S;

/* Вычисляем интеграл с половинным шагом */
h2 = h; h /= 2.0;
if( h == 0.0 ) break; /* потеря значимости */
n = n2; n2 *= 2;

Seven = Seven + Sodd;
/* Вычисляем сумму по новым точкам: */
for( Sodd = 0.0, x = a+h, i = 0;
i < n;
i++, x += h2 )
Sodd += (*f)(x);

/* Значение интеграла */
S = h / 3.0 * (fab + 4.0 * Sodd + 2.0 * Seven );

} while( niter < 31 && fabs(S - Sprev) / 15.0 >= eps );
/* Используем условие Рунге для окончания итераций */

return ( 16.0 * S - Sprev ) / 15.0 ;
/* Возвращаем уточненное по Ричардсону значение */
}