Maximum area enclosure given side lengths

A peer of mine gave me the following problem :

Given a sequence of $n$ lengths (i.e.,$L_1, L_2, .., L_n$ ) where each is the length of the side, find a sequence of $n$ points (where $p_k = (x_k, y_k)$) such that $dist(p_k, p_{k+1}) = L_k$ and $dist(p_1, p_n) = L_n$ where $dist(a, b)$ is the Euclidean distance between points $a$ and $b$.

Note that the first point , $p1$, must be $(0, 0)$ and the second point , $p_2$, must be $(0, L_1)$ .

These points must correspond to the ordered clockwise vertices of the simple polygon having the maximum possible area for the given side lengths.

Determine the coordinates of the points that correspond to a polygon of maximal area.

Now, I've seen some similar problems, but this one I don't even know how to begin with. Determining the maximum area is quite easily done by Brahmagupta's formula, but the actual finding of the points seems really hard.

The only thing I could think of was finding some sort of enclosing circle and then trying to determine the points, but that's not even a starting point.

Any ideas on how to solve this?

• Your starting point is correct: the maximum area is obtained by a cyclic polygon (all vertices lie on a common circle). I don't know a way to find the radius of the circle except to form a nonlinear equation to solve for it in terms of subtended angles that sum to $2\pi$. Sep 17, 2016 at 18:47
• Yup, finding the radius and center seem to be the actual problem. Sep 17, 2016 at 18:53
• This is a litteral copy-paste od a question on a rated online contest that is still in progress. Why are you cheatting? Sep 18, 2016 at 17:21
• @plamenko: There is a site-policy to lock contest questions under certain conditions. Can you point out where this problem is posed online, or give more information? Sep 18, 2016 at 18:22
• @hardmath Here is the link hackerrank.com/contests/w23/challenges/enclosure Sep 18, 2016 at 19:11

(As far as I can see, the online competition is now over, so I shall reinstate my answer as promised.)

As commented by hardmath, and agreed by the OP, the solution is a cyclic polygon whose vertices lie on a common circle.

Note that if one of the line segments would be longer than the sum of the length of the rest of the line segments, they cannot form a polygon. Because of this, we know that $$L_i \le \sum_{k \ne i}^{n} L_k$$ for all $i = 1 \ldots n$. (This also means that $0 \le \theta_i \lt \pi$.)

Each of the $n$ line segments essentially forms an isosceles triangle with the center of the common circle. Two of the sides are of length $r$, the radius of the common circle, and the third is the line segment itself, length $L_i$. (Essentially, we can decompose the polygon into triangular wedges, each being an isosceles triangle. This also means that the order of the line segments does not affect the area, as reordering the triangular parts does not change their areas.)

The length of a circular chord fulfills $$L_i = 2 r \sin\left(\frac{\theta_i}{2}\right)$$ Here, $\theta_i$ is the angle between the $r$-length edges in the isosceles triangles. Solving for $\theta_i$, we see that $$\theta_i = 2 \arcsin \left( \frac{L_i}{2 r} \right) = \arccos \left( 1 - \frac{ L_i^2 }{2 r^2} \right)$$ because $2 \arcsin(x) = \arccos(1 - 2 x^2)$ when $0 \le x \le 1$.

Note that this also means that $$\frac{d \, \theta_i}{d \, r} = - \frac{2 L}{\sqrt{4 r^4 - L^2 r^2}}$$

In order for the figure to be a polygon, the sum of the angles $\theta_i$ must equal $180°$: $$\sum_{i=1}^{n} \theta_i = 2 \pi$$

Combining the above two, and dividing by two, we get $$\sum_{i=1}^{n} \arcsin \left( \frac{L_i}{2 r} \right) = \pi$$ which we can use to numerically solve for $r$.

In order for each $\theta_i$ to exist, $\max(L)/(2r) \le 1$. This means that $r \ge \max(L)/2$, and that $r_{min} = \max(L)/2$ is the minimum allowed radius in the numerical methods.

We also know that the radius must be less than half the perimeter, so $r_{max} = (\sum_i L_i) / 2$.

When we have found the $r$ for which $\sum_i \theta_i = 2 \pi$, we can compute the vertex coordinates. Since the points are to be in clockwise order, we can use $$\begin{cases} x_i' = x_0 - r \cos \varphi_i \\ y_i' = y_0 + r \sin \varphi_i \\ \end{cases}$$ where $$\begin{cases} x_0 = r \cos \varphi_1 \\ y_0 = -r \sin \varphi_1 \end{cases}$$ and $$\varphi_i = \begin{cases} -\frac{\theta_i}{2}, & i = 1 \\ -\frac{\theta_i}{2} + \sum_{k=1}^{i-1} \theta_i, & 2 \le i \le n \end{cases}$$ The above gives us $\varphi_1 = -\theta_1/2$ and $\varphi_2 = -\varphi_1 = \theta_1/2$, so the first line segment will be vertical; $x_1 = x_2 = y_1 = 0$, $y_2 \gt 0$ (in fact $y_2 = L_1$).

The above therefore yields $x_1 = 0$, $y_1 = 0$, $x_2 = 0$, $y_2 = L_1$, with the rest of the vertices $(x_i, y_i)$ in clockwise order around the center at $(x_0, y_0)$.

Here is an example program that reads line lengths from standard input, and uses a binary search to find the radius with at least six digits of precision, outputting the vertex coordinates to standard output. example.c:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define  PI     3.141592653589793238462643383279502884197
#define  TWO_PI 6.283185307179586476925286766559005768394

double theta(const double L, const double r)
{
return 2.0 * asin(0.5 * L / r);
}

double d_theta(const double L, const double r)
{
const double r2 = r * r;
return -2.0 * L / (r2 * sqrt(4.0 - L * L / r2));
}

double find_radius(const double length[], const size_t count, const double epsilon, size_t *iteration_count)
{
const double min_theta = 1.0 - epsilon;
const double max_theta = 1.0 + epsilon;

double sum_length = length[0];
double max_length = length[0];
size_t iterations = 0;
size_t i;

for (i = 1; i < count; i++) {
sum_length += length[i];
if (max_length < length[i])
max_length = length[i];
}

while (1) {
double sum_theta = 0.0;

iterations++;

for (i = 0; i < count; i++)

sum_theta /= TWO_PI;

if (sum_theta >= min_theta && sum_theta <= max_theta)
break;
else
if (sum_theta < 1.0)
else
}

if (iteration_count)
*iteration_count = iterations;

}

int main(void)
{
double  *line_length   = NULL;
size_t   line_count    = 0;
size_t   line_maxcount = 0;
double   radius, phi, x, y, x0, y0;
size_t   i;

/* Read line lengths from standard input. */
while (1) {
char   buffer[1024], *line;
double length;

line = fgets(buffer, sizeof buffer, stdin);
if (!line)
break;

if (sscanf(line, " %lf", &length) < 1)
break;

if (length <= 0.0)
break;

if (line_count >= line_maxcount) {
const size_t maxcount = (line_count | 255) + 257;
double      *temp;

temp = realloc(line_length, maxcount * sizeof line_length[0]);
if (!temp) {
fprintf(stderr, "Out of memory.\n");
return EXIT_FAILURE;
}

line_length = temp;
line_maxcount = maxcount;
}

line_length[line_count++] = length;
}

/* Verify sanity. */
if (line_count >= 3) {
double max_length = line_length[0];
double sum_length = line_length[0];

for (i = 1; i < line_count; i++) {
sum_length += line_length[i];
if (max_length < line_length[i])
max_length = line_length[i];
}

if (max_length > sum_length - max_length) {
fprintf(stderr, "Not a valid polygon; one of the line segments is too long.\n");
return EXIT_FAILURE;
}
} else {
fprintf(stderr, "Need at least three line lengths.\n");
return EXIT_FAILURE;
}

fprintf(stderr, "Read %zu line lengths:\n", line_count);
for (i = 0; i < line_count; i++)
fprintf(stderr, "\t%.6f\n", line_length[i]);

phi = -0.5 * theta(line_length[0], radius);

for (i = 0; i < line_count; i++) {
x = x0 - radius * cos(phi);
y = y0 + radius * sin(phi);
printf("%.6f %.6f\n", x, y);
}

return EXIT_SUCCESS;
}


Compile it using e.g. in Linux

gcc -Wall -O2 example.c -lm -o example


and run it, generating say 10 random line lengths as input, using

awk 'BEGIN { srand(); for (i=0; i<10; i++) printf "%.6f\n", 10.0*rand() }' | ./example


Here is an example of the output, first 12 lines to standard error, followed by 10 lines to standard output:

Read 10 line lengths:
1.487793
5.241667
2.002344
4.271945
7.815320
5.452837
6.012328
7.506031
3.599992
1.636063
radius = 7.377844 using 19 iterations.
0.000000 0.000000
0.000000 1.487793
2.346549 6.174880
3.990793 7.317614
8.195604 8.071990
14.300092 3.191984
14.080471 -2.256429
9.609570 -6.276271
2.286102 -4.630880
0.344424 -1.599406

• I tried newtons method to solve for r but the convergence is too slow. Are their any faster methods? Sep 18, 2016 at 3:07
• The speed of Newton's method depends on the adequacy of intial guesses. Once you are in the "basin of attraction" of the root (of this univariate problem, which lacks root multiplicity), convergence should be rapid (quadratic convergence). Sep 18, 2016 at 3:30
• @hardmath I am using $max(L_k)/2+c$ as the starting point where $max(L_k)$ enables the $arcsin$ to be defined at all places and c is a small constant like .01 Sep 18, 2016 at 4:10
• @AmanDeepGautam: If you put some sample lengths in a Comment on my Answer, I'll show how a more robust root-finding method can narrow the estimates enough to use Newton's method effectively. Sep 18, 2016 at 4:16
• I've removed my answer pending a resolution of the comment on the Question that it "cut-and-pastes" an ongoing contest problem. Sep 18, 2016 at 18:55

I'll add some remarks about the setup and solution for circumscribing diameter $d$, from which the maximum area (of the cyclic polygon) can be obtained.

(0) The case $n=3$ is directly solved by using Heron's formula for area of the triangle. The case $n=4$ also has a direct solution, but it is less well-known. For $n \gt 4$ we present an iterative approach (solution of a nonlinear equation).

(1) The order of the chord (polygon side) lengths does not affect the maximum area that the polygon can attain. There is a nice proof of this by transposition of adjacent chords, showing that any permutation of sides in a cyclic polygon preserves area.

(2) Even if one specifies an order of the sides lengths $L_1,L_2,\ldots,L_n$, we can without loss of generality arrange (by rotation) that $L_1$ is the maximum of these. Then the check $L_1 \lt \sum_{i=2}^n L_i$ guarantees the existence of a polygon with positive area.

(3) The diameter $d$ of the circle must be at least the maximum chord length $L_1$. There are two distinguishable cases where the diameter strictly exceeds $L_1$. Either the vertices all lie on a semicircle (and the circle's center is outside the polygon), or else the center of the circle will lie in the interior of the polygon.

(4) A good first step to solving for the radius (equiv. the diameter) of the circle distinguishes which of those arrangements is feasible. For the sides $L_i$ with $i\gt 1$, define the sum of half the angles they will subtend in a circle of (yet undetermined) diameter $d$:

$$\Theta(d) := \sum_{i=2}^n \sin^{-1} (L_i/d)$$

(5) If $\Theta(L_1) = \pi/2$, then $d=L_1$ is the solution, and besides the diameter chord $L_1$, the other chords are inscribed in its semicircle.

(6) If $\Theta(L_1) \lt \pi/2$, then the feasible $d$ will exceed $L_1$ but all the vertices will lie on a semicircle, with the longest chord $L_1$ bridging across the smaller chords. The nonlinear equation to satisfy is then:

$$\Theta(d) = \sin^{-1} (L_1/d)$$

(7) If $\Theta(L_1) \gt \pi/2$, then the feasible $d$ will exceed $L_1$ but not all the vertices will lie on a semicircle, and thus the circle's center will lie in the interior of the cyclic polygon. Thus central angles subtend by the chords will thus add to $2\pi$, and in terms of our $\Theta(d)$ notation, the nonlinear equation to solve is:

$$\Theta(d) = \pi - \sin^{-1} (L_1/d)$$

Example

Consider the lengths, arranged as above so that the longest length comes first:

$$10, 8, 2, 5, 6, 4, 8, 5, 4, 2$$

The check $\Theta(10) \approx 4.771 \gt \pi/2$ tells us the vertices will wind around the center of the circle, so that per (7) the nonlinear equation to satisfy is:

$$\sin^{-1}(L_1/d) + \Theta(d) = \pi$$

It simplifies the derivatives (needed for Newton's method) if we make the change of variable $z = L_1/d$, so the equation becomes:

$$F(z) := \sin^{-1}(z) + \sum_{i=2}^n \sin^{-1}(\alpha_i z) = \pi$$

where $\alpha_i = L_i/L_1$. For $z\in (0,1)$ we get its derivative:

$$F'(z) = (1-z^2)^{-1/2} + \sum_{i=2}^n \alpha_i(1 - \alpha_i^2 z^2)^{-1/2}$$

Note that $F(z)$ is continuous, monotone increasing on $[0,1]$. Now $F(0) = 0$, and we also know as a result of the check $\Theta(10) \gt \pi/2$ that $F(1) \gt \pi$. Therefore a root of $F(z) = \pi$ will exist "in between" by IVT.

Often one uses a step or two of a low-order convergent method such as bisection or secant method to "get close" to the root before applying a Newton method. A bisection step amounts to checking $F(1/2) \approx 2.760 \lt \pi$, so our root search should be restricted to $z\in [1/2,1]$.

Since $F(1/2) \approx 2.760$ and $F(1) \approx 6.342$, a step of the secant method (linear interpolation) gives a root approximation near:

$$z_0 = \frac{(1/2)(F(1)-\pi) - (1)(F(1/2)-\pi)}{F(1)-F(1/2)} \approx 0.5533$$

which we adopt as an initial guess for Newton's method:

$$z_{k+1} = z_k - \frac{F(z_k) - \pi}{F'(z_k)}$$

A table of a few subsequent iterates, calculated using Google Sheets, shows quick convergence:

$$\begin{array}{c|c|c|c} k & z_k & F(z_k) & F'(z_k) \\ \hline 1 & 0.5533006452 & 3.070530159 & 5.884332671 \\ 2 & 0.5653772049 & 3.141749930 & 5.910562711 \\ 3 & 0.5653505955 & 3.141592654 & 5.910503695 \\ 4 & 0.5653505954 & 3.141592654 & 5.910503695 \end{array}$$

For Newton's method it is typical that iterates quickly settles down to approaching the limit from one side (as they did in this case).

Recall that $z = L_1/d$ and thus $d = L_1/z \approx 17.68813915$. Converting from diameter to radius gives $r \approx 8.844069575$.

• The following edges list $2$ $4$ $5$ $8$ $4$ $6$ $5$ $2$ $8$ $10$ takes around 13000 iterations to converge for a precision of $10^{-5}$. I take 5.1 as the start value. Could you please tell how to take better start points. Sep 18, 2016 at 4:26
• The other thing I realized was that with these values, the estimated radius oscillates around the actual value. So for the above case the actual value is around $8.8$ then the method quickly (in around 10 iterations) realizes the value in range $[8.1, 9.3]$, but then it takes forever to reduce that interval range. For example the values of estimated radius starting iteration 12 from my program for above example are $9.36149, 8.26006, 9.3421, 8.28434, 9.32485$ and so on Sep 18, 2016 at 12:48
• @AmanDeepGautam: The oscillating behavior is a clue that something is wrong in the coding. Newton's method will rarely jump back and forth around the limit, but rather iterates will tend to the limit from one consistent side in cases where the second derivative is of one sign around the root. Sep 18, 2016 at 18:18
• Very good points, hardmath. In my answer, I tried to show (1) and (2) by decomposing the polygon into wedges made of isosceles triangles, with two sides of $r$, and the third $L_i$. The rule used to fix $r$ ensures the wedges form a complete polygon. It is obvious that their order is then irrelevant, as $r$ only scales the wedges/triangles. I like your explanation and points much better than mine, however. I haven't examined the parameter space for $r$ yet to see how smooth or not it is, and whether there are many local minima and maxima. I think I'll write a test program to find out. Sep 18, 2016 at 18:18