$ontext
# Maximum area for unit-diameter polygon (Prieto's second model).

# This derives from a GAMS model by Francisco J. Prieto -- the second
# one in E-mail that he sent to Margaret Wright.  Prieto described
# this as "a very particular solution for the problem, taking advantage
# of [a] certain property that I [Prieto] am convinced (although I cannot
# prove it) that the solution must have.  The formulation
# is much more efficient, and you can go up to sizes of around 70-80.

# "Symmetry is assumed (so the size is reduced to a half), and also as
# many constraints as possible are set up as equalities.

# "The results obtained seem to coincide (within some error margin) with
# the ones obtained from the previous program."
$offtext

$if not setglobal N $setglobal N 26
abort$(%N% gt 2*floor(%N%/2)) "N must be even";

set vertices /1*%N%/;
alias(j,k,vertices);

scalar Lsup;
Lsup = %N%/2 - 1;

set i(vertices);
i(vertices) = yes$(ord(vertices) le Lsup);

variables area;
positive variables x(j), y(j);
x.up(j) = 1; y.up(j) = 1;
x.l(j) = .5;
y.l(vertices)$i(vertices) = ord(vertices)/Lsup;

equation obj, ordered(k), diam(k,j), extr;

* distance constraints (tight)
diam(k,j)$(i(k) and ord(k) ge Lsup/2 and ord(j) ge Lsup-ord(k) and ord(j) le ord(k) and ord(j) le Lsup+1-ord(k))..
        sqr(x(k) + x(j)) + sqr(y(k) - y(j)) =e= 1;

* extra distance constraint
extr.. sum(j$(ord(j) eq Lsup), sqr(x(j)) + sqr(y(j))) =e= 1;

* second coordinate constraint (ordered values)
ordered(k)$(i(k) and ord(k) gt 1)..
  y(k) =g= y(k-1);

obj..
  area =e= sum(j$(ord(j) lt Lsup), y(j)*(x(j-1)-x(j+1))) +
                sum(j$(ord(j) eq Lsup), x(j) + y(j)*x(j-1));

model pgon /obj,ordered,diam,extr/;
solve pgon using nlp max area;

option area:8;
display area.l;
display x.l, y.l;