Approximate Square Formation
Pyth, 28 bytes
jbc*\#Qho.a,,N*NJ_/_QN,JQ_SQ

Try it online.
The crux is that I sort potential m on the following property:
(m - ceil(k/m))^2 + (m*ceil(k/m) - k)^2

Note the total absence of n. The total shape is defined merely by m. Then I transform the above property once more, and my final sorting weight is defined as the Euclidean distance between the following two points:
(m, m*ceil(k/m)) and (ceil(k/m), k)

This changes the weight values, but not their ordering.