Here’s the “proper solution”:
// Problem: integrate a two dimensional function "f(x,y)" using limits:
// "y" from 0 to 1 and "x" from "y" to 1.
// Solution: change the variable "x:[y, 1]" -> "t:[0, 1]" and
// then the integration goes from 0 to 1 in both dimensions.
// So "x = (1 - y) * t + y" and "dx/dt = (1 - y)".
// See also:
// http://en.wikipedia.org/wiki/Integration_by_substitution
// http://en.wikipedia.org/wiki/Multiple_integral#Change_of_variables
//
double f2(const double *X)
{
double f = 0.0;
double t = X[0]; // "t" (i.e. the original "x")
double y = X[1]; // "y"
double x = (1.0 - y) * t + y; // calculate new "x"
// calculate the original function "f(x,y)"
f = x * y; // f(x,y) = x * y
// f = x * y * y; // f(x,y) = x * y^2
f *= (1.0 - y); // multiply the result by "dx/dt"
return f;
}
Just for reference: “Non-rectangular region integration”