#include <kv/doubleint-singular.hpp>

struct Func {
	template <class T> T operator() (const T& x, const T& y) {
		return sqrt(x * cos(y)) * (cos(x * y));
	}
};

struct Func_f {
	template <class T> T operator() (const T& x, const T& y) {
		return x * cos(y);
	}
};

struct Func_g {
	template <class T> T operator() (const T& x, const T& y) {
		return cos(x * y);
	}
};

typedef kv::interval<double> itv;

int main()
{
	std::cout.precision(17);

	std::cout << kv::doubleintegral_power3(Func_f(), Func_g(), itv(0.), itv(0.125), itv(0.), itv(0.125), 12, itv(0.5), 1, 0) << "\n";
}