test/math/simpson.cpp


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

#include "../util.h"
std::function<double(double)> f;
constexpr double EPS = 1e-9;
#include <math/simpson.cpp>

struct RandomPolynom {
	vector<double> polynom;
	RandomPolynom(int deg) : polynom(deg) {
		for (auto& x : polynom) x = Random::real<double>(-100, 100);
	}
	double operator()(double x) const {
		double res = 0;
		double xx = 1;
		for (double y : polynom ) {
			res += xx * y;
			xx *= x;
		}
		return res;
	}
	double area(double a, double b) const {
		double res = 0;
		double aa = a;
		double bb = b;
		ll d = 1;
		for (double y : polynom) {
			res += bb / d * y;
			res -= aa / d * y;
			aa *= a;
			bb *= b;
			d++;
		}
		return res;
	}
};

void stress_test() {
	timer t;
	ll queries = 0;
	for (int tries = 0; tries < 1'000; tries++) {
		ll n = Random::integer<ll>(0, 6);
		RandomPolynom poly(n);
		f = poly;
		for (ll i = 0; i < 200; i++) {
			double l = Random::real<double>(-20, 20);
			double r = Random::real<double>(-20, 20);
			if (l > r) swap(l, r);

			t.start();
			double got = integrate(l, r);
			t.stop();
			double expected = poly.area(l, r);
			if (float_error(got, expected) > 1e-6) cerr << fixed << setprecision(20) << "got: " << got << ", expected: " << expected << FAIL;
			queries++;
		}
	}
	if (t.time > 5000) cerr << "too slow: " << t.time << FAIL;
	cerr << "tested random queries: " << queries << " (" << t.time << "ms)" << endl;
}

int main() {
	stress_test();
}