#include <stdio.h>
#include <math.h>
// from http://r...content-available-to-author-only...e.org/wiki/Numerical_integration#C
double int_leftrect(double from, double to, double n, double (*func)())
{
double h = (to-from)/n;
double sum = 0.0, x;
for(x=from; x <= (to-h); x += h)
sum += func(x);
return h*sum;
}
double integral(double (*func)()) {
return int_leftrect(0, 1, 10, func);
}
static double from;
static double to;
static double (*fn)();
double scaled(double x) {
return fn(from + x * (to - from));
}
double integral_range(double (*func)(), double a, double b) {
from = a;
to = b;
fn = func;
return integral(scaled) * (b - a);
}
double test_fn(double x) {
double result = 2 * x;
printf("TRACE: test_fn(%f) => %f\n", x
, result
); return result;
}
int main(void) {
double result = integral_range(test_fn, 3, 9);
double expected = (9 * 9) - (3 * 3);
printf("result of numerical integration: %f\n", result
); printf("error of numerical integration: %f\n", fabs(result
- expected
)); return 0;
}
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgovLyBmcm9tIGh0dHA6Ly9yLi4uY29udGVudC1hdmFpbGFibGUtdG8tYXV0aG9yLW9ubHkuLi5lLm9yZy93aWtpL051bWVyaWNhbF9pbnRlZ3JhdGlvbiNDCmRvdWJsZSBpbnRfbGVmdHJlY3QoZG91YmxlIGZyb20sIGRvdWJsZSB0bywgZG91YmxlIG4sIGRvdWJsZSAoKmZ1bmMpKCkpCnsKICAgZG91YmxlIGggPSAodG8tZnJvbSkvbjsKICAgZG91YmxlIHN1bSA9IDAuMCwgeDsKICAgZm9yKHg9ZnJvbTsgeCA8PSAodG8taCk7IHggKz0gaCkKICAgICAgc3VtICs9IGZ1bmMoeCk7CiAgIHJldHVybiBoKnN1bTsKfQoKZG91YmxlIGludGVncmFsKGRvdWJsZSAoKmZ1bmMpKCkpIHsKCXJldHVybiBpbnRfbGVmdHJlY3QoMCwgMSwgMTAsIGZ1bmMpOwp9CgoKc3RhdGljIGRvdWJsZSBmcm9tOwpzdGF0aWMgZG91YmxlIHRvOwpzdGF0aWMgZG91YmxlICgqZm4pKCk7CmRvdWJsZSBzY2FsZWQoZG91YmxlIHgpIHsKCXJldHVybiBmbihmcm9tICsgeCAqICh0byAtIGZyb20pKTsKfQoKZG91YmxlIGludGVncmFsX3JhbmdlKGRvdWJsZSAoKmZ1bmMpKCksIGRvdWJsZSBhLCBkb3VibGUgYikgewoJZnJvbSA9IGE7Cgl0byA9IGI7CglmbiA9IGZ1bmM7CglyZXR1cm4gaW50ZWdyYWwoc2NhbGVkKSAqIChiIC0gYSk7Cn0KCmRvdWJsZSB0ZXN0X2ZuKGRvdWJsZSB4KSB7Cglkb3VibGUgcmVzdWx0ID0gMiAqIHg7CglwcmludGYoIlRSQUNFOiB0ZXN0X2ZuKCVmKSA9PiAlZlxuIiwgeCwgcmVzdWx0KTsKCXJldHVybiByZXN1bHQ7Cn0KCmludCBtYWluKHZvaWQpIHsKCWRvdWJsZSByZXN1bHQgPSBpbnRlZ3JhbF9yYW5nZSh0ZXN0X2ZuLCAzLCA5KTsKCWRvdWJsZSBleHBlY3RlZCA9ICg5ICogOSkgLSAoMyAqIDMpOwoJcHJpbnRmKCJyZXN1bHQgb2YgbnVtZXJpY2FsIGludGVncmF0aW9uOiAlZlxuIiwgcmVzdWx0KTsKCXByaW50ZigiZXJyb3Igb2YgbnVtZXJpY2FsIGludGVncmF0aW9uOiAlZlxuIiwgZmFicyhyZXN1bHQgLSBleHBlY3RlZCkpOwoJcmV0dXJuIDA7Cn0K