Approximate Square Formation Background

I have a bunch of square-shaped boxes of equal size, and since I'm a neat person, I want to arrange them all into a square formation.
However, their number is not necessarily a perfect square, so I may have to approximate the square shape.
I want you to find me the most aesthetically pleasing arrangement -- programmatically, of course.

Input

Your input is a single positive integer k, representing the number of boxes.

Output

Your program shall choose two positive integers m, n such that m*(n-1) &lt; k &le; m*n holds.
They represent the width and height of the large square-like shape we are arranging.
Since we are looking for aestethically pleasing shapes, the quantity (m - n)2 + (m*n - k)2 shall be minimal, so that the shape is close to a square, and its area is close to k.
If there are still several candidates for the pair (m, n), choose the one where the width m is maximal.

Now, your actual output shall not be the numbers m and n.
Instead, you shall print the arrangement of boxes, using the character # to represent a box.
More specifically, you shall print n-1 rows, each of which consists of m characters #, and then one row of k - m*(n-1) characters #.
Note that the output contains exactly k characters #.

Rules and Scoring

There shall not be any leading or trailing whitespace in the output, except that the last row may be padded with trailing spaces to be of length m, if desired.
There may be one trailing newline, but no preceding newlines.
You may use any printable ASCII character in place of #, if desired.

You may write a full program, or return a string from a function.
The lowest byte count wins, and standard loopholes are disallowed.

Test Cases

Here are the correct outputs for a few input values.

1
#
2
##
3
##
#
4
##
##
8
###
###
##
13
#####
#####
###
17
######
######
#####
18
#####
#####
#####
###
20
#####
#####
#####
#####
21
######
######
######
###
22
######
######
######
####
23
#####
#####
#####
#####
###
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.