11
\$\begingroup\$

Given consecutive side lengths s1, s2, s3... s_n of an n-gon inscribed in a circle, find its area. You may assume that the polygon exists. In addition, the polygon will be convex and not self-intersecting, which is enough to guarantee uniqueness. Built-ins that specifically solve this challenge, as well as built-in functions that calculate the circumradius or circumcenter, are banned (this is different from a previous version of this challenge).

Input: The side lengths of the cyclic polygon; may be taken as parameters to a function, stdin, etc.

Output: The area of the polygon.

The answer should be accurate to 6 decimal places and must run within 20 seconds on a reasonable laptop.

This is code golf so shortest code wins!

Specific test cases:

[3, 4, 5] --> 6
[3, 4, 6] --> 5.332682251925386
[3, 4, 6, 7] --> 22.44994432064365
[5, 5, 5, 5] --> 25
[6, 6, 6, 6, 6] --> 61.93718642120281
[6.974973020933265, 2.2393294197257387, 5.158285083300981, 1.4845682771595603, 3.5957940796134173] --> 21.958390804292847
[7.353566082457831, 12.271766915518073, 8.453884922273897, 9.879017670784675, 9.493366404245332, 1.2050010402321778] --> 162.27641678140589

Test case generator:

function randPolygon(n) {
  var left = 2 * Math.PI;
  var angles = [];
  for (var i = 0; i < n - 1; ++i) {
    var r = Math.random() * left;
    angles.push(r);
    left -= r;
  }
  angles.push(left);
  var area = 0;
  var radius = 1 + Math.random() * 9;
  for (var i = 0; i < angles.length; ++i) area += radius * radius * Math.sin(angles[i]) / 2;
  var sideLens = angles.map(function(a) {
    return Math.sin(a / 2) * radius * 2;
  });

  document.querySelector("#radius").innerHTML = radius;
  document.querySelector("#angles").innerHTML = "[" + angles.join(", ") + "]";
  document.querySelector("#inp").innerHTML = "[" + sideLens.join(", ") + "]";
  document.querySelector("#out").innerHTML = area;

  draw(angles);
}

function draw(angles) {
  var canv = document.querySelector("#diagram"),
    ctx = canv.getContext("2d");
  var size = canv.width
  ctx.clearRect(0, 0, size, size);
  ctx.beginPath();
  ctx.arc(size / 2, size / 2, size / 2, 0, 2 * Math.PI, true);
  ctx.stroke();
  ctx.beginPath();
  ctx.moveTo(size, size / 2);
  var runningTotal = 0;
  for (var i = 0; i < angles.length; ++i) {
    runningTotal += angles[i];
    var x = Math.cos(runningTotal) * size / 2 + size / 2;
    var y = Math.sin(runningTotal) * size / 2 + size / 2;
    ctx.lineTo(x, y);
  }
  ctx.stroke();
}
document.querySelector("#gen").onclick = function() {
  randPolygon(parseInt(document.querySelector("#sideLens").value, 10));
}
<div id="hints">
  <p><strong>These are to help you; they are not part of the input or output.</strong>
  </p>
  Circumradius:
  <pre id="radius"></pre>
  Angles, in radians, of each sector (this are NOT the angles of the polygon):
  <pre id="angles"></pre>
</div>
<hr>
<div id="output">
  Input:
  <pre id="inp"></pre>
  Output:
  <pre id="out"></pre>
</div>
<hr>
<div id="draw">
  Diagram:
  <br />
  <canvas id="diagram" width="200" height="200" style="border:1px solid black"></canvas>
</div>

Number of side lengths:
<input type="number" id="sideLens" step="1" min="3" value="3" />
<br />
<button id="gen">Generate test case</button>

\$\endgroup\$
4
  • 7
    \$\begingroup\$ I know an easy way to find its perimeter. \$\endgroup\$
    – mIllIbyte
    Apr 16, 2016 at 18:04
  • 1
    \$\begingroup\$ I know an easy way to find the number of sides \$\endgroup\$
    – Luis Mendo
    Apr 16, 2016 at 19:05
  • \$\begingroup\$ This problem is pretty easy given the circumradius, but without it it's incredibly difficult. \$\endgroup\$
    – poi830
    Apr 16, 2016 at 19:32
  • \$\begingroup\$ It's also easy if there are fewer than five sides, not that it matters in code golf. \$\endgroup\$
    – Neil
    Apr 16, 2016 at 19:54

3 Answers 3

5
\$\begingroup\$

Python 2, 195 191 bytes

from math import*
C=sorted(input());l,h=0,1/max(C)
while h-l>2e-16:m=(l+h)/2;a=[asin(c*m)for c in C[:-1]];f=pi-sum(a);l,h=[l,m,h][sin(f)/m>C[-1]:][:2]
print sum(sin(2*t)/l/l for t in a+[f])/8

Uses a binary search to find the inverse diameter, then calculates the area of each segment by the angle/radius.

It finds the radius by first summing all but the largest chord angle, and checking the remaining angle to the remaining chord. Those angles are then also used to compute the area of each segment. A segment's area can be negative, if it's angle is bigger than 180 degrees. We do a binary search over the inverse diameter so the domain is bounded and well-behaved.

Readable implementation:

def segment_angles(line_segments, invd):
    return [2*math.asin(c*invd) for c in line_segments]

def cyclic_ngon_area(line_segments):
    line_segments = list(sorted(line_segments))
    lo, hi = 0, 1/max(line_segments)
    while hi - lo > 2e-16:
        mid = (lo + hi) / 2
        angles = segment_angles(line_segments[:-1], mid)
        angles.append(2*math.pi - sum(angles))
        if math.sin(angles[-1]/2) / mid > line_segments[-1]:
            lo = mid
        else:
            hi = mid
    return sum([math.sin(a)/lo/lo/8 for a in angles])
\$\endgroup\$
7
  • \$\begingroup\$ Does this work if the center is outside the polygon? (For example, a triangle with side lengths 6, 7, 12). Sometimes the sqrt(4**2 - c**2/4) needs to be negative, when the angle is greater than pi. \$\endgroup\$
    – soktinpk
    Apr 16, 2016 at 19:58
  • \$\begingroup\$ @soktinpk I fixed my answer. \$\endgroup\$
    – orlp
    Apr 20, 2016 at 8:43
  • \$\begingroup\$ The code gives wrong results for the input [1000000, 1000000, 1999999]. The exact result for this test is 1999999*sqrt(3999999)/4 = 999999375.0000546875... \$\endgroup\$ Oct 25, 2022 at 19:03
  • \$\begingroup\$ @HugoPfoertner I fixed it by squaring the initial hi estimate, at the cost of 3 bytes. \$\endgroup\$
    – orlp
    Oct 26, 2022 at 0:46
  • 1
    \$\begingroup\$ @HugoPfoertner Right, it's unbounded. I fixed my answer by instead doing a binary search over the inverse diameter, which lies in (0, 1/max(edges)]. There's still numerical precision loss, but it should give 6 significant digits of precision for reasonable-sized inputs. \$\endgroup\$
    – orlp
    Oct 26, 2022 at 11:57
5
\$\begingroup\$

PARI/GP, 182 bytes

A(v)=n=#s=vecsort(v);g=1^s;h=sum(k=1,n-1,asin(s[k]/S=s[n]))<Pi/2;g[n]-=2*h;d=solve(D=S,S/(vecsum(s)/2/S-1),sum(k=1,n,g[k]*asin(s[k]/D))-Pi*!h);sum(k=1,n,g[k]*s[k]*sqrt(d^2-s[k]^2))/4

Using @alephalpha's hints that global variables can be used and his improved chained assignments, significantly shorter code is possible. Also an unnecessary factor 8 in calculation of Dmax was removed.

Version with comments

AreaCyclicPolygon(VectorofSidelenghts) =
{
  my( s = vecsort(VectorofSidelenghts), \\ allows access to longest side
                                        \\ and loop over shorter sides
      n = #s,
      S = s[n],    \\ assignment postponed into loop in golfed code
      g = vector(n,i,1), \\ Auxiliary vector of ones
      d,
      h = sum ( k=1, n-1, asin(s[k]/S)) < Pi/2, \\ determine whether the triangle 
                                                   \\ with longest side makes an angle > Pi.
      Dmax = S/sqrt(vecsum(s)/2/S - 1) \\ upper limit of diameter from asymptotic
                                             \\ expansion showing 1/sqrt(eps) dependency.
                                             \\ In "golfed" version sqrt can be ommitted
                                             \\ drastically overestimating Dmax; solve still works
    );
  g[n] -= 2*h; \\ factor = -1 for largest triangle
               \\ when the polygon only occupies less than half of the circle
  d = solve ( D=s[n], Dmax,
              sum (k=1, n, g[k]*asin(s[k]/D)) - Pi*!h); \\ solve sum of angles equation
                                                        \\ including special case of
                                                        \\ large angle for longest side
  sum (k=1, n, g[k]*s[k]*sqrt(d^2 - s[k]^2))/4  \\ return value: use solution for d in sum
                                                \\ of triangle areas, again with potentially
                                                \\ negative contribution of triangle
                                                \\ with longest side
};

Examples

AreaCyclicPolygon([3,4,5])
 6.000000000000000000
AreaCyclicPolygon([3,4,6])
 5.3326822519253854377215978803358864040
AreaCyclicPolygon([0.25,0.25,0.499999])
 0.00012499968750010937502343751074219391159
AreaCyclicPolygon([5,5,5,5])
 25.000000000000000000000000000000000000
AreaCyclicPolygon([6,6,6,6,6])
 61.937186421202809219324431185989945578
AreaCyclicPolygon([1e6,1e6,2e6-1])
 999999375.00005468750292968783569307190

Note on estimating the maximum diameter

Regardless of whether the diameter (or equivalently the radius) is calculated using bisection or with a function for solving a non-linear equation (like PARI's solve here), a bracketing interval is required for the solution value. Trivially, the diameter cannot be smaller than the longest side. However, since the polygon can degenerate when the sum of the shorter sides approaches the longest side, the diameter of the polygon can become arbitrarily large. In order to estimate this diameter, it is sufficient to use a triangle as a model. From the well-known Heron formula for the area A one gets the radius $$ r(a,b,c) = \frac{a b c}{\sqrt{(a+b-c) (a-b+c) (-a+b+c) (a+b+c)}} $$ If one assumes a model triangle \$\left(\frac{1}{4},\frac{1}{4},\frac{1}{2}-\epsilon \right)\$ normalized to an approximate perimeter of 1 with a scaled length deficit \$\epsilon\$ of the longest side, then r becomes $$ r = \frac{\frac{1}{2}-\epsilon }{16 \sqrt{\left(\frac{1}{2}-\epsilon \right)^2 (1-\epsilon ) \epsilon }} $$ Assuming \$\epsilon\$ to be small, r can be expanded r in a series $$ r = \frac{1}{16 \sqrt{\epsilon }}+\frac{\sqrt{\epsilon }}{32}+\frac{3 \epsilon ^{3/2}}{128}+O\left(\epsilon ^{5/2}\right) $$ The leading term contains the singularity \$\frac{1}{\sqrt{\epsilon }}\$, which must be taken into account when calculating the upper bound. It doesn't matter whether you use the circumference or the longest side to scale back.

Whether one computes \$\epsilon\$ as $$ \epsilon = \frac{a+b+c}{2 c}-1 = \frac{a+b-c}{2 c} $$ which is the vecsum(s)/2/S - 1 in the program, or as $$ \epsilon = 1-\frac{2 c}{a+b+c} $$ makes no essential difference for small \$\epsilon\$.

\$\endgroup\$
5
  • \$\begingroup\$ Welcome to Code Golf, and nice answer! \$\endgroup\$ Oct 23, 2022 at 22:17
  • 1
    \$\begingroup\$ Welcome to Code Golf! I'm really happy to see a new PARI/GP answer. Using global variables is fine in code golf. So you don't need to use my. my(s=vecsort(v),n=#v,S=s[n],g=vector(n,i,1),d,h=sum(k=1,n-1,asin(s[k]/S))<Pi/2) can become S=s[n=#s=vecsort(v)];g=1^s;h=sum(k=1,n-1,asin(s[k]/S))<Pi/2. \$\endgroup\$
    – alephalpha
    Oct 24, 2022 at 11:55
  • \$\begingroup\$ Thank you for the hints. I don't know enough about the rules here. May I incorporate your suggestions into my code? As an editor in OEIS, I constantly look at PARI code and occasionally create code myself. But there we often have the opposite problem, that we criticize poorly readable code. But it's fun here and I plan to contribute from time to time. \$\endgroup\$ Oct 24, 2022 at 12:13
  • \$\begingroup\$ Since my reputation is not yet high enough to make comments on other posts, I'll say this at this point: @orlp 's Python code returns an incorrect result with the input of my last test example, because the starting value for the upper bound "hi" is set much too low for this input. print (cyclic_ngon_area ([1000000, 1000000, 1999999])) -> 124019623218.91309. Setting the hi value correctly is more laborious and requires longer code. \$\endgroup\$ Oct 25, 2022 at 8:16
  • \$\begingroup\$ @HugoPfoertner Sorry. When I said you can use global variables, I meant that you can use global variables in functions (instead of local variables declared by my), not that you can assume the input is stored in a global variable. You should still define a function that takes the input. \$\endgroup\$
    – alephalpha
    Oct 27, 2022 at 23:48
0
\$\begingroup\$

Octave, 89 bytes

r=sum(s=input(''));while sum(a=asin(s/2/r))<pi r*=1-1e-4;b=a;end;disp(sum(cos(b).*s/2*r))

Explanation

The angle a spanned by a segment of length s is 2*asin(s/2/r), given a circumradius r. Its area is cos(a)*s/2*r.

Algorithm

  1. Set r to something too large, such as the perimeter.
  2. If the cumulated angle is smaller than 2pi, reduce r and repeat step 2.
  3. Calculate the area.
\$\endgroup\$
10
  • \$\begingroup\$ On average, how many iterations does this take for r to be set? (out of curiosity) \$\endgroup\$
    – soktinpk
    Apr 16, 2016 at 19:50
  • \$\begingroup\$ There is no way this has the required precision. You repeatedly multiply the radius by 0.9999 to reduce it, this makes it very easy to undershoot the required 6 decimals of precision. \$\endgroup\$
    – orlp
    Apr 16, 2016 at 19:53
  • \$\begingroup\$ @soktinpk around 15000 for r*=1-1e-4 and 150000 for r*=1-1e-5. \$\endgroup\$
    – Rainer P.
    Apr 16, 2016 at 19:54
  • \$\begingroup\$ @RainerP. Those two values are the same. \$\endgroup\$
    – anon
    Apr 16, 2016 at 19:55
  • 1
    \$\begingroup\$ @soktinpk it is generally not a good idea to make an exception for a specific answer. \$\endgroup\$
    – Cyoce
    Apr 16, 2016 at 20:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.