I thought I knew how to do this. Took all day.

permute_formula_variables.cpp

//maps permutations of decreasing ranges (EG, the first number may be total, the second number can only be up to total-1, the third only up to total-2, etc ) to non-overlapping numbers (none of the numbers in tri will be the same after this has been run. Used for combinatorial generation. Spaced must only contain falses (or nothing). It must be large enough for the largest allowable possible value in tris
void triangularToNonoverlapping(vector<bool>& spaced, vector<uint>& tris, uint total){
	for(uint i=0; i<tris.size(); ++i){
		uint place = tris[i];
		uint nextPlace = place;
		uint scani = 0;
		while(scani <= place){
			if(spaced[scani]){
				place += 1;
			}
			scani += 1;
		}
		spaced[place] = true;
		tris[i] = place;
	}
}

bool iterateTriangular(vector<uint>& tris, uint total){ //false iff past end. If that happens, don't iterate again. It cycles.
	assert(tris.size() <= total);
	uint i=tris.size()-1;
	again:
		if(tris[i] >= total - 1 - i){
			tris[i] = 0;
			if(i == 0){
				return false;
			}else{
				--i;
				goto again;
			}
		}
		tris[i] += 1;
	return true;
}

//iterates the vars.size() - ownedStart last variables over every disjoint combination of the points in dimensions
template<typename F>
void overPermutationsOfCoords(Coord dimensions, vector<VarIdentification>& vars, uint ownedStart, F f){
	uint total = dimensions.x*dimensions.y;
	uint dif = vars.size() - ownedStart;
	uint remainingTotal = total - dif;
	vector<bool> space;
	vector<uint> triangledValueAssignments;
	vector<uint> valueAssignments;
	triangledValueAssignments.resize(dif, 0);
	valueAssignments.resize(dif, 0);
	while(true){
		space.clear();
		space.resize(total, false);
		//write positions of previous values to space
		for(uint i=0; i<ownedStart; ++i){
			space[vars[i].valueNumber] = true;
		}
		//mingle the new ones into space, and make them nonoverlapping instead of triangular (I was surprised to find this could be done at the same time, but the algorithm was so simple)
		triangularToNonoverlapping(space, valueAssignments, remainingTotal);
		//update the values
		for(uint i=0; i<valueAssignments.size(); ++i){
			VarIdentification& vi = vars[ownedStart + i];
			vi.valueNumber = valueAssignments[i];
			int iy = vi.valueNumber/dimensions.x;
			vi.value = Coord(vi.valueNumber - iy*dimensions.x, iy);
		}
		
		//yeild
		bool brb = false;
		f(brb);
		if(brb){ return; }
		
		//iterate the permutation
		if(!iterateTriangular(triangledValueAssignments, remainingTotal)){ return; }
		valueAssignments = triangledValueAssignments;
	}
}