#include <kv/doubleint-singular.hpp>

struct Func {
	template <class T> T operator() (const T& x, const T& y) {
		return pow((1-cos(x)) * cos(y), 1/T(3)) * (cos(x * y));
	}
};

struct Func_f {
	template <class T> T operator() (const T& x, const T& y) {
		return (1-cos(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, 1./itv(3.), 2, 0) << "\n";
}