The smallest area of a convex grid polygon

Question

I got an email from Hugo Pfoertner, an Editor-in-Chief at the On-Line Encyclopedia of Integer Sequences, with a terrific idea for a fastest-code challenge, which will also help verify or expand the boundary of an open mathematical question in optimization.

Here's the idea: we want to confirm and extend OEIS sequences A063984 and A070911, which are about finding convex polygons on the grid \$\mathbb Z \times \mathbb Z\$ that have \$n\$ vertices and minimal area. (See an example below.) These sequences had been computed up to \$n=25\$, but unfortunately, some of the programs that were used to compute these sequences have been lost, and it would be useful to have new programs that we could use to confirm the sequence.

If you are able to extend the sequence or correct any erroneous terms, Hugo or I will update the OEIS to give you credit for that—plus, you'll get lots of brownie points from me. :)

Challenge

Your program should take a positive integer n (\$n \geq 3\$) and output a set of \$n\$ vertices of an area-minimal convex \$n\$-gon with vertices in the square grid. It should also output twice¹ the area of the corresponding polygon, which should be equal to (or less than!) \$A070911(n)\$.

This is a fastest-code challenge, and the score will be based on the largest \$n\$ that you can compute on my machine² in under \$10\$ minutes. In case of a tie, whichever program computes the next value fastest will win.

It would be helpful for me if you would post the runtimes and results from running your program on your machine, so that I can get an order-of-magnitude estimate for how quickly your program runs.

Examples

The OEIS sequence A070911 begins (starting at the index \$n=3\$): $$1, 2, 5, 6, 13, 14, 21, 28, 43, 48, 65, 80, 103, 118, 151, 174, 213, 242, 289, 328, 387, 420, 497, \dots$$

Hugo Pfoertner posted a PDF on the OEIS containing the best known examples. I have included the example for \$n=13\$ here, which shows that the minimum (known) area for a convex lattice \$13\$-gon is \$\displaystyle \frac{65}{2}\$.

If you want to play around with some examples, Markus Sigg wrote a visualizer in Javascript (unfortunately, no HTTPS for this one).

---

¹ Why twice the area? So that all of the solutions are whole numbers.

² I'll be running this on my 2021 MacBook Pro on an Apple M1 Pro chip and with 16 GB of RAM.

Answer 1

JavaScript (Node.js), \$n=19\$

This is a very early attempt. Given 10 minutes, it is only able to compute up to \$n=19\$ on my laptop and finds the same values as the ones listed in A070911.

It assumes \$n\ge7\$ in order to avoid some edge cases. Although the code could be updated to support \$n<7\$, I'm not sure it's worth the effort.

A few assumptions are made about the shape of the \$n\$-gon and the size of its enclosing rectangle. The later one should probably be relaxed to make sure that no shorter solution is missed.

Each solution comes with an ascii-art representation.

Code

const MAX_N    = 13;
const MAX_TIME = 600;

let pathCache = {},
    ts = Date.now(),
    min = 4, max = 4;

for(let n = 7; n <= MAX_N; n++) {
  console.log(`n = ${n}`);

  let [ score, solution ] = solve(n, min, max),
      w = solution.w - 1,
      h = solution.h - 1;

  console.log(`Area * 2 = ${score}`);
  console.log(`Enclosing rectangle: ${w} x ${h}`);
  displaySolution(solution);

  let time = Math.round((Date.now() - ts) / 1000);

  console.log(`Total time: ${time}s\n`);

  if(time >= MAX_TIME) {
    break;
  }
  min = solution.h;
  max = solution.w + 2;
}

function solve(n, min, max) {
  let part = partitionsInFour(n),
      rect = [],
      best = Infinity,
      solution;

  for(let w = min; w <= max; w++) {
    for(let h = min; h <= w; h++) {
      rect.push([ w, h ]);
    }
  }

  for(let [ w, h ] of rect) {
    for(let y0 = 0; y0 < h - 1; y0++) {
      for(let x0 = 0; x0 < w - 1; x0++) {
        for(let y1 = 1; y1 < h; y1++) {
          for(let x1 = 1; x1 < w; x1++) {
            let [ score, res ] = solveRectangle(part, w, h, y0, x0, y1, x1);

            if(score < best) {
              best = score;
              solution = { w: w, h: h, y0: y0, res: res };
            }
          }
        }
      }
    }
  }
  return [ best, solution ]
}

function solveRectangle(part, w, h, y0, x0, y1, x1) {
  let best = Infinity,
      solution;

  for(let p of part) {
    let p0, p1, p2, p3;

    if(
      (p0 = getPaths(x0,         h - y0 - 2)[p[0] - 1]) &&
      (p1 = getPaths(w - x0 - 2, h - y1 - 1)[p[1] - 1]) &&
      (p2 = getPaths(w - x1 - 1, y1 - 1    )[p[2] - 1]) &&
      (p3 = getPaths(x1 - 1,     y0        )[p[3] - 1])
    ) {
      let outerArea = p0.area + p1.area + p2.area + p3.area,
          score = (w - 1) * (h - 1) * 2 - outerArea;

      if(score < best) {
        best = score;
        solution = [
          p0.path,
          [...p1.path].reverse(),
          p2.path,
          [...p3.path].reverse()
        ];
      }
    }
  }
  return [ best, solution ];
}

function getPaths(w, h) {
  if(pathCache[[ w, h ]]) {
    return pathCache[[ w, h ]];
  }

  let set = new Set(), list = [];

  function search(x = 0, y = 0, pw = 0, ph = 1, area = 0, path = []) {
    if(x == w && y == h) {
      let len = path.length,
          outerArea = w * h * 2 - area;

      if(!set.has(len) || outerArea > list[len].area) {
        set.add(len);

        list[len] = {
          area: outerArea,
          path: path
        }
      }
    }

    for(let w0 = 1; w0 <= w - x; w0++) {
      for(let h0 = 1; h0 <= h - y; h0++) {
        if(h0 * pw < w0 * ph) {
          search(
            x + w0, y + h0, w0, h0,
            area + ((w - x) * 2 - w0) * h0, [ ...path, [ w0, h0 ]]
          );
        }
      }
    }
  }

  search();
  pathCache[[ w, h ]] = list;
  pathCache[[ h, w ]] = list.map(o => ({
    area: o.area,
    path: [...o.path].reverse().map(([ dx, dy ]) => [ dy, dx ])
  }));

  return list;
}

function partitionsInFour(n) {
  let list = [];

  (function search(n, i = 1, l = []) {
    if(n) {
      if(i <= n && l.length != 4) {
        search(n - i, 1, [i, ...l]);
        search(n, i + 1, l);
      }
    }
    else if(l.length == 4) {
      list.push(l);
    }
  })(n);

  return list;
}

function displaySolution(solution) {
  let m = [...Array(solution.h)].map(_ => [...Array(solution.w)].fill('--')),
      x, y, i = 1;

  function mark(x, y) {
    if(m[y][x] == '--') {
      m[y][x] = i.toString().padStart(2, '0');
    }
    i++;
  }

  mark(x = 0, y = solution.y0);
  mark(x, ++y);
  solution.res[0].forEach(([ dx, dy ]) => mark(x += dx, y += dy));
  mark(++x, y);
  solution.res[1].forEach(([ dx, dy ]) => mark(x += dx, y -= dy));
  mark(x, --y);
  solution.res[2].forEach(([ dx, dy ]) => mark(x -= dx, y -= dy));
  mark(--x, y);
  solution.res[3].forEach(([ dx, dy ]) => mark(x -= dx, y += dy));

  console.log(m.map(r => r.join(' ')).join('\n'));
}

Try it online! (up to \$n=13\$)

Output

n = 7
Area * 2 = 13
Enclosing rectangle: 3 x 3
01 07 -- --
02 -- -- 06
-- -- -- 05
-- 03 04 --
Total time: 0s

n = 8
Area * 2 = 14
Enclosing rectangle: 3 x 3
-- 08 07 --
01 -- -- 06
02 -- -- 05
-- 03 04 --
Total time: 0s

n = 9
Area * 2 = 21
Enclosing rectangle: 4 x 4
01 09 -- -- --
02 -- -- 08 --
-- -- -- -- 07
-- 03 -- -- 06
-- -- 04 05 --
Total time: 0s

n = 10
Area * 2 = 28
Enclosing rectangle: 5 x 4
-- 10 09 -- -- --
01 -- -- -- 08 --
02 -- -- -- -- 07
-- 03 -- -- -- 06
-- -- -- 04 05 --
Total time: 0s

n = 11
Area * 2 = 43
Enclosing rectangle: 6 x 5
-- 11 10 -- -- -- --
01 -- -- -- -- 09 --
02 -- -- -- -- -- 08
-- -- -- -- -- -- 07
-- 03 -- -- -- 06 --
-- -- 04 05 -- -- --
Total time: 1s

n = 12
Area * 2 = 48
Enclosing rectangle: 6 x 6
-- 12 11 -- -- -- --
01 -- -- -- 10 -- --
02 -- -- -- -- 09 --
-- -- -- -- -- -- --
-- 03 -- -- -- -- 08
-- -- 04 -- -- -- 07
-- -- -- -- 05 06 --
Total time: 4s

n = 13
Area * 2 = 65
Enclosing rectangle: 8 x 6
-- 13 12 -- -- -- -- -- --
01 -- -- -- -- 11 -- -- --
02 -- -- -- -- -- -- 10 --
-- -- -- -- -- -- -- -- 09
-- 03 -- -- -- -- -- -- 08
-- -- 04 -- -- -- -- 07 --
-- -- -- -- 05 06 -- -- --
Total time: 7s

n = 14
Area * 2 = 80
Enclosing rectangle: 8 x 7
-- -- -- 13 12 -- -- -- --
-- 14 -- -- -- -- 11 -- --
01 -- -- -- -- -- -- 10 --
02 -- -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- 09
-- 03 -- -- -- -- -- -- 08
-- -- 04 -- -- -- -- 07 --
-- -- -- -- 05 06 -- -- --
Total time: 22s

n = 15
Area * 2 = 103
Enclosing rectangle: 9 x 8
-- -- 14 13 -- -- -- -- -- --
-- 15 -- -- -- -- 12 -- -- --
-- -- -- -- -- -- -- -- 11 --
01 -- -- -- -- -- -- -- -- 10
02 -- -- -- -- -- -- -- -- 09
-- -- -- -- -- -- -- -- -- --
-- 03 -- -- -- -- -- -- 08 --
-- -- 04 -- -- -- -- 07 -- --
-- -- -- -- 05 06 -- -- -- --
Total time: 38s

n = 16
Area * 2 = 118
Enclosing rectangle: 9 x 9
-- -- -- -- 14 13 -- -- -- --
-- -- 15 -- -- -- -- 12 -- --
-- 16 -- -- -- -- -- -- 11 --
-- -- -- -- -- -- -- -- -- --
01 -- -- -- -- -- -- -- -- 10
02 -- -- -- -- -- -- -- -- 09
-- -- -- -- -- -- -- -- -- --
-- 03 -- -- -- -- -- -- 08 --
-- -- 04 -- -- -- -- 07 -- --
-- -- -- -- 05 06 -- -- -- --
Total time: 70s

n = 17
Area * 2 = 151
Enclosing rectangle: 11 x 10
-- -- -- 16 15 -- -- -- -- -- -- --
-- 17 -- -- -- -- -- 14 -- -- -- --
01 -- -- -- -- -- -- -- -- 13 -- --
02 -- -- -- -- -- -- -- -- -- 12 --
-- -- -- -- -- -- -- -- -- -- -- --
-- 03 -- -- -- -- -- -- -- -- -- 11
-- -- -- -- -- -- -- -- -- -- -- 10
-- -- -- -- -- -- -- -- -- -- -- --
-- -- -- 04 -- -- -- -- -- -- 09 --
-- -- -- -- 05 -- -- -- -- 08 -- --
-- -- -- -- -- -- 06 07 -- -- -- --
Total time: 98s

n = 18
Area * 2 = 174
Enclosing rectangle: 12 x 10
-- -- -- -- 16 15 -- -- -- -- -- -- --
-- -- 17 -- -- -- -- -- 14 -- -- -- --
-- 18 -- -- -- -- -- -- -- -- 13 -- --
-- -- -- -- -- -- -- -- -- -- -- 12 --
01 -- -- -- -- -- -- -- -- -- -- -- --
02 -- -- -- -- -- -- -- -- -- -- -- 11
-- -- -- -- -- -- -- -- -- -- -- -- 10
-- 03 -- -- -- -- -- -- -- -- -- -- --
-- -- 04 -- -- -- -- -- -- -- -- 09 --
-- -- -- -- 05 -- -- -- -- -- 08 -- --
-- -- -- -- -- -- -- 06 07 -- -- -- --
Total time: 202s

n = 19
Area * 2 = 213
Enclosing rectangle: 13 x 12
-- -- 18 17 -- -- -- -- -- -- -- -- -- --
-- 19 -- -- -- -- 16 -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- 15 -- -- -- -- --
01 -- -- -- -- -- -- -- -- -- -- -- -- --
02 -- -- -- -- -- -- -- -- -- -- 14 -- --
-- -- -- -- -- -- -- -- -- -- -- -- 13 --
-- -- -- -- -- -- -- -- -- -- -- -- -- --
-- 03 -- -- -- -- -- -- -- -- -- -- -- 12
-- -- -- -- -- -- -- -- -- -- -- -- -- 11
-- -- 04 -- -- -- -- -- -- -- -- -- -- --
-- -- -- 05 -- -- -- -- -- -- -- -- 10 --
-- -- -- -- -- 06 -- -- -- -- -- 09 -- --
-- -- -- -- -- -- -- -- 07 08 -- -- -- --
Total time: 436s

Answer 2

Fortran n=17, 323 s

C Find strictly convex polygons of minimum area on square grid
C Author: Hugo Pfoertner, 2018
      implicit integer (A-Z)
      parameter (vlarge=2147483647)
C Number of vertices
      parameter (n=17)
C Number of coordinate pairs to be used as polygon edges,
C read from external list 
      parameter (ms=2048)
C Length of lists for segments, coordinates, areas
C Must be extended for n>nm
      parameter (nm=25)
      dimension xd(ms), yd(ms), x(nm), y(nm), a(nm), as(nm), nn(nm)
      equivalence
     &  (n1,nn(1)),(n2,nn(2)),(n3,nn(3)),(n4,nn(4)),(n5,nn(5)),
     &  (n6,nn(6)),(n7,nn(7)),(n8,nn(8)),(n9,nn(9)),(n10,nn(10)),
     &  (n11,nn(11)),(n12,nn(12)),(n13,nn(13)),(n14,nn(14)),
     &  (n15,nn(15)),(n16,nn(16)),(n17,nn(17)),(n18,nn(18)),
     &  (n19,nn(19)),(n20,nn(20)),(n21,nn(21)),(n22,nn(22)),
     &  (n23,nn(23)),(n24,nn(24)),(n25,nn(25))
C Number of polygons with minimal area found
      integer*8 count
C File names of external files, command argument
      character*15 bspirx, bspiry, carg
C Progress indicator line
      character pline*150
C CPU time
      real cptime
C function to calculate d^2 of enclosing circle,
C to be replaced by function encirc
C if exact enclosing circle is needed
      integer diamet
      external diamet
C variables  needed in diameter calculation
      doubleprecision xc, yc, rc, d, diamin, diamax
C Some choices for OEIS files describing spirals
C Square spiral
C      data bspirx, bspiry /'b174344.txt', 'b274923.txt' /
C Circular rings (Sloane)
C      data bspirx, bspiry /'b283307.txt', 'b283308.txt' /
C Circular rings
      data bspirx, bspiry /'../b305575.txt', '../b305576.txt' /
C Statement function: Double area of triangle
      triar(x1,y1, x2,y2, x3,y3) =
     &      x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2)
C Progress indicator
      pline = '....+....1....+....2....+....3....+....4....+....5' //
     &  '....+....6....+....7....+....8....+....9....+....A' //
     &  '....+....B....+....C....+....D....+....E....+....F'
C
C read external files with coordinates of points in spiral
C X
      open ( unit=10, file=bspirx, status='old',
     &       form='formatted', iostat=ios)
      if ( ios .ne. 0 ) stop 'Error opening bfile spiral x'
      do 1 i = 1, ms
      read (10,*,end=999) k, xd(i)
1     continue
      close (unit=10)
C Y
      open ( unit=10, file=bspiry, status='old',
     &       form='formatted', iostat=ios)
      if ( ios .ne. 0 ) stop 'Error opening bfile spiral y'
      do 2 i = 1, ms
      read (10,*,end=999) k, yd(i)
2     continue
      close (unit=10)
C
C For convenience: write first nonnegative (x,y) pairs to terminal
      do 3 i = 2, 120
      if (xd(i) .ge. 0 .and. yd(i) .ge. 0 ) write (*,1003)i,xd(i),yd(i)
1003  format ( 3 i3 )
3     continue
C
C preset minimum area
      ami = vlarge
C if an upper bound is known: least area + 1
C      ami = 183
C preset diameter extreme values
      diamin = 1.0D20
      diamax = 0.0D0
C Total number of polygons with same minimum area
      count = 0
C get number of list items from first parameter of program call
      CALL get_command_argument(1, carg)
      read (carg, *) m
      if ( m .gt. ms ) stop 'm exceeds length of segment list'
      write (*,*) 'Segments used:', m
C get index of
      CALL get_command_argument(2, carg)
      read (carg, *) n2first
      if ( xd(n2first) .lt. 0  .or. yd(n2first) .lt. 0 )
     &   stop 'illegal negative start step'
C limit range, assuming first coordinate pair on files is (0,0)
      n2first = max(2,min (m,n2first))
      write (*,*) 'First start step:', xd(n2first), yd(n2first)
C
C Start building the polygon
C
C Freeze first point
      x(1) = 0
      y(1) = 0
      n1 = 0
C
C loop over second point
      do 20 n2 = n2first, m
      L = 2
      x(L) = xd(nn(L))
      y(L) = yd(nn(L))
C
C Limit to angle 0 <= Pi/2
      if ( x(L) .lt. 0 .or. y(L) .lt. 0 .or. y(L) .gt. x(L) ) goto 20
C
C optional: Exclusion of extremely long segments
C      if (dble(x(L)**2 + y(L)**2) .gt. diamin) goto 20
      write (*,1006) count, n2, xd(nn(L)), yd(nn(L))
1006  format (/,'Min area polygons found so far: ', i0,
     & ', next n2 = ', i0, ' (',i0,',',i0,')')
C if wanted: progress indicator 
      WRITE(*, 1004, ADVANCE='NO') pline(1:1)
1004  format (A1)

C loop over third point
      do 30 n3 = 2, m
C Progress indicator
      WRITE(*, 1004, ADVANCE='NO') pline(n3:n3)
      L = 3
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      as(L) = a(L)
      if ( a(L) .le. 0 ) goto 30
      if ( a(L) .gt. ami-n+L ) goto 30
C The following blocks are repeated in code with adaptation for
C current segment number (code easily generated by a small script
C or preprocessor)
      do 40 n4 = 2, m
      L = 4
C try extension by segment from list
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
C area contribution
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 40
C left turn?
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0) goto 40
C start point still left of straight line through endpoints of segment?
      if ( triar(x(1),y(1), x(2),y(2), x(L),y(L)) .le. 0 ) goto 40
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 40
C
      do 50 n5 = 2, m
      L = 5
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar ( x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 50
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0) goto 50
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 50
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 50
C
      do 60 n6 = 2, m
      L = 6
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 60
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0) goto 60
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 60
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 60
C
      do 70 n7 = 2, m
      L = 7
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 70
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0) goto 70
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 70
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 70
C
      do 80 n8 = 2, m
      L = 8
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 80
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0) goto 80
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 80
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 80
C
      do 90 n9 = 2, m
      L = 9
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 90
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0) goto 90
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 90
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 90
C
      do 100 n10 = 2, m
      L = 10
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 100
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0)goto 100
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 100
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 100
C 
      do 110 n11 = 2, m
      L = 11
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 110
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0)goto 110
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 110
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 110
C
      do 120 n12 = 2, m
      L = 12
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 120
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0)goto 120
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 120
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 120
C
      do 130 n13 = 2, m
      L = 13
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 130
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0)goto 130
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 130
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 130
C
      do 140 n14 = 2, m
      L = 14
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 140
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0)goto 140
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 140
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 140
C
      do 150 n15 = 2, m
      L = 15
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 150
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0)goto 150
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 150
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 150
C
      do 160 n16 = 2, m
      L = 16
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 160
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0)goto 160
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 160
      as(L) = as(L-1) + a(L)
      if ( as(L) .gt. ami-n+L ) goto 160
C
      do 170 n17 = 2, m
      L = 17
      x(L) = x(L-1) + xd(nn(L))
      y(L) = y(L-1) + yd(nn(L))
      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
      if ( a(L) .le. 0 ) goto 170
      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0)goto 170
      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 170
      as(L) = as(L-1) + a(L)
C in last line of repeated code part n=L
c      if ( as(L) .gt. ami-n+L ) goto 170
C example how to continue for n>17
c      do 180 n18 = 2, m
c      L = 18
c      x(L) = x(L-1) + xd(nn(L))
c      y(L) = y(L-1) + yd(nn(L))
c      a(L) = triar (x(1),y(1), x(L-1),y(L-1), x(L),y(L))
c      if ( a(L) .le. 0 ) goto 180
c      if (triar(x(L-2),y(L-2), x(L-1),y(L-1), x(L),y(L)) .le. 0)goto 180
c      if ( triar (x(L),y(L), x(1),y(1), x(2),y(2)) .le. 0 ) goto 180
c      as(L) = as(L-1) + a(L)
c      if ( as(L) .gt. ami-n+L ) goto 180
C ...
C ...
C
C Update minimum
      if ( as(L) .lt. ami ) then
        count = 0
        ami = as(L)
C
C alternative with exact enclosing circle
C        call encirc ( 1, n, x, y, xc, yc, rc )
C        diamin = 4*rc**2
        call cpu_time (cptime)
C type cast assumed to work diamin (real*8) = diamet (integer)
        diamin = diamet (n,x,y)
        write (*,1000) n, as(L), diamin, (x(k),y(k),k=1,n)
1000    format (/,i2, 1X, i0, f14.6, 2x, 25('(',i0,',',i0,')',:,',') )
        write (*,1001) cptime, nn(2:n)
1001    format ( F12.4,' s: ',i0, 25(1X,i0,:) )
      endif
C
C check for multiple solutions with same mimimum area
      if ( as(L) .eq. ami ) then
        call cpu_time (cptime)
        d = diamet(n,x,y)
C        call encirc ( 1, n, x, y, xc, yc, rc )
C        d = 4*rc**2
        count = count + 1
        if ( d .lt. diamin ) then
          diamin = d
          write (*,1000) n, as(L), diamin, (x(k),y(k),k=1,n)
          write (*,1001) cptime, nn(2:n)
        endif
        if ( d .gt. diamax ) then
          diamax = d
          write (*,1000) n, as(L), -diamax, (x(k),y(k),k=1,n)
          write (*,1001) cptime, nn(2:n)
        endif
      endif
250   continue
240   continue
230   continue
220   continue
210   continue
200   continue
190   continue
180   continue
170   continue
160   continue
150   continue
140   continue
130   continue
120   continue
110   continue
100   continue
90    continue
80    continue
70    continue
60    continue
50    continue
40    continue
30    continue
20    continue
C
      call cpu_time ( cptime )
      write (*,1007) cptime, count
1007  format (/,'CPU time: ', f12.4, ' s',/,
     &        'Number of polygons with minimum area: ', i0)
999   continue
      end
C
C Maximum of mutual point distance sufficient as an estimate.
C Exact enclosing circle needs a more sophisticated method,
C e.g., Welz's algorithm
      integer function diamet (n, x, y)
      integer n, x(*), y(*)
      id = 0
      do 10 i = 1, n-1
      do 20 j = i+1, n
      jd = (x(i)-x(j))**2 + (y(i)-y(j))**2
      id = max (id,jd)
20    continue
10    continue
      diamet = id
      end

This is essentially the first version of the program I started with in 2018 just for illustration with no tweaks. It's more to have a place to make some general notes on pitfalls of this problem that I think are important. When those are scattered in comments on individual answers, it's hard to keep track of. The program only handles the case n=17, which I chose because it is the smallest n without a proof of optimality. In order to run the program, 2 external files are required, namely the b-files of the OEIS sequences A305575 and A305576, which should be one directory above the program. A typical output looks like this:

.\17.exe 56 1
  2  1  0
  3  0  1
...
114  6  1
115  1  6
 Segments used:          56
 First start step:           1           0

Min area polygons found so far: 0, next n2 = 2 (1,0)
...
17 185    373.000000  (0,0),(1,0),(1,1),(0,3),(-1,4),(-4,6),(-6,7),(-9,8),(-13,9),(-14,9),(-16,8),(-17,7),(-17,6),(-16,5),(-13,3),(-11,2),(-8,1)
      0.0469 s: 2 3 16 7 41 17 33 53 4 18 8 5 9 45 21 37

17 185   -373.000000  (0,0),(1,0),(1,1),(0,3),(-1,4),(-4,6),(-6,7),(-9,8),(-13,9),(-14,9),(-16,8),(-17,7),(-17,6),(-16,5),(-13,3),(-11,2),(-8,1)
      0.0469 s: 2 3 16 7 41 17 33 53 4 18 8 5 9 45 21 37

17 178    370.000000  (0,0),(1,0),(1,1),(0,3),(-2,5),(-5,7),(-7,8),(-10,9),(-14,10),(-15,10),(-16,9),(-16,8),(-15,6),(-14,5),(-11,3),(-9,2),(-6,1)
      0.0938 s: 2 3 16 23 41 17 33 53 4 8 5 20 9 45 21 37

17 159    306.000000  (0,0),(1,0),(1,1),(0,3),(-2,6),(-3,7),(-5,8),(-8,9),(-12,10),(-13,10),(-14,9),(-14,8),(-13,6),(-12,5),(-9,3),(-7,2),(-4,1)
      0.1094 s: 2 3 16 40 7 17 33 53 4 8 5 20 9 45 21 37

17 159    296.000000  (0,0),(1,0),(1,1),(0,3),(-2,6),(-3,7),(-6,9),(-8,10),(-11,11),(-12,11),(-13,10),(-13,9),(-12,7),(-10,4),(-9,3),(-7,2),(-4,1)
      0.1094 s: 2 3 16 40 7 41 17 33 4 8 5 20 44 9 21 37

17 157    265.000000  (0,0),(1,0),(1,1),(0,4),(-1,6),(-3,9),(-4,10),(-6,11),(-9,12),(-10,12),(-11,11),(-11,10),(-10,7),(-9,5),(-8,4),(-5,2),(-3,1)
      0.2969 s: 2 3 32 16 40 7 17 33 4 8 5 36 20 9 45 21
.+.
17 157    232.000000  (0,0),(1,0),(2,1),(2,2),(1,4),(0,5),(-3,7),(-5,8),(-8,9),(-9,9),(-11,8),(-12,7),(-12,6),(-11,4),(-10,3),(-8,2),(-5,1)
      0.9844 s: 2 6 3 16 7 41 17 33 4 18 8 5 20 9 21 37

17 151    202.000000  (0,0),(1,0),(2,1),(2,2),(1,4),(-1,7),(-2,8),(-5,10),(-7,11),(-8,11),(-9,10),(-10,8),(-10,7),(-9,5),(-8,4),(-5,2),(-3,1)
      1.0469 s: 2 6 3 16 40 7 41 17 4 8 19 5 20 9 45 21

17 151    193.000000  (0,0),(1,0),(2,1),(2,2),(1,5),(0,7),(-1,8),(-3,9),(-6,10),(-7,10),(-9,9),(-10,8),(-10,7),(-9,5),(-8,4),(-5,2),(-3,1)
      1.1406 s: 2 6 3 32 16 7 17 33 4 18 8 5 20 9 45 21

17 151    137.000000  (0,0),(1,0),(2,1),(3,3),(3,4),(2,7),(1,9),(0,10),(-2,11),(-3,11),(-5,10),(-6,9),(-7,7),(-7,6),(-6,4),(-5,3),(-2,1)
      2.1562 s: 2 6 15 3 32 16 7 17 4 18 8 19 5 20 9 45
...1....+....2....+....3....+....4....+....5....+.
Min area polygons found so far: 24, next n2 = 6 (1,1)
....+....1....+....2....+....3....+....4....+....5....+.
Min area polygons found so far: 48, next n2 = 10 (2,0)
....+....1....+....2....+....3....+....4....+....5....+.
Min area polygons found so far: 48, next n2 = 14 (2,1)
....+....1....+....2....+....3....+....4....+....5....+.
Min area polygons found so far: 65, next n2 = 22 (2,2)
....+....1....+....2....+....3....+....4....+....5....+.
Min area polygons found so far: 65, next n2 = 26 (3,0)
....+....1....+....2....+....3....+....4....+....5....+.
Min area polygons found so far: 65, next n2 = 30 (3,1)
....+....1....+....2....+....3....+....4....+....5....+.
Min area polygons found so far: 70, next n2 = 38 (3,2)
....+....1....+....2....+....3....+....4....+....5....+.
Min area polygons found so far: 76, next n2 = 46 (4,0)
....+....1....+....2....+....3....+....4....+....5....+.
Min area polygons found so far: 76, next n2 = 50 (4,1)
....+....1....+....2....+....3....+....4....+....5....+.
CPU time:     323.0781 s
Number of polygons with minimum area: 76

For n=17 there is not only the compact solution found by everyone, but also exotic needle-shaped solutions, like Squared diameter 1361. This is the one I've known so far with the largest diameter. As far as I know, it has not been proven that no extreme solutions of this kind with a smaller area exist. If one could show that there are no other solutions with even greater stretching, then that would be an important step towards a proof of optimality for this n.

Heavily stretched polygons

I don't want to spoil anyone's good mood, but if you all only find the solutions that I gave 4 years ago, the doubts remain whether we are making things too easy for ourselves when searching. I have already pointed out the existence of very strongly stretched polygons with likewise minimal areas. As a test case for your programs, you can try to find at least one even slimmer solution than the following at n=19: n=19, 2*Area=213, Diameter^2=1105 My program finds the shown and -intentionally undisclosed- slimmer solutions (squared diameters = 5*prime number, prime, ..) after about 200 s (17700 s for the prime squared diameter) using 1000 points from the spiral files.

If other programs also find these solutions, then that would increase my confidence considerably.

Update

In the meantime I have found that all of the strongly distorted and needle-shaped polygons found by my programs can be reduced to the already known slightly deformed shapes by applying shear transformations. So far I haven't found any exceptions to this observation. Apparently, allowing longer polygon sides does not bring any advantage in terms of further reducing the area. In a way, this contradicts the asymptotic elliptical shape described in the work of Bárány, I., Tokushige, N. (2004) or version without paywall with semi-axes \$a=0.003573 n^2\$ and \$b=1.656 n\$, which at \$n=27\$ gives an axis ratio of about \$1/15\$. The observed ratio of the solution with \$2*A(27)=625\$, which has meanwhile also been found by my own program, is only about \$1/2\$.

Fortran revised, n=17, 8.8 s

See Lattice-Polygons on GitHub. Another faster version exists (approx. 70% of run time), with explicit expansion of the inner loops, similar to code shown above, but the GitHub version is cleaner.

Answer 3

Python, n=44 in under 2 minutes

This can go up to 44-gons in under 2 minutes, using pypy on my laptop. CPython is slightly slower (2'45'').

The approach is rigorous. It uses the dynamic programming algorithm in the following paper:

• David Eppstein, Mark Overmars, Günter Rote und Gerhard Woeginger, Finding minimum area k-gons. Discrete and Computational Geometry 7 (1992), 45-58. https://doi.org/10.1007/BF02187823

• see also: Joseph S. B. Mitchell, Günter Rote, Gopalakrishnan Sundaram, and Gerhard Woeginger, Counting convex polygons in planar point sets. Information Processing Letters 56 (1995), 45-49. http://page.mi.fu-berlin.de/rote/Papers/abstract/Counting+k-subsets+and+convex+k-gons+in+the+plane.html

The algorithm is also explained in the comments of the code.

Here are some smallest polygons:

• 20-gon, area=242/2: [(0,0),(3,1),(5,2),(6,3),(7,5),(7,6),(6,8),(5,9),(3,10),(0,11),(-1,11),(-4,10),(-6,9),(-7,8),(-8,6),(-8,5),(-7,3),(-6,2),(-4,1),(-1,0)]
• 21-gon, area=289/2: [(0,0),(2,1),(5,3),(6,4),(8,7),(9,9),(10,12),(10,13),(9,14),(7,15),(6,15),(3,14),(1,13),(-2,11),(-3,10),(-4,8),(-5,5),(-5,4),(-4,2),(-3,1),(-1,0)]
• 22-gon, area=328/2: [(0,0),(3,1),(5,2),(8,4),(9,5),(10,7),(10,8),(9,10),(8,11),(6,12),(3,13),(2,13),(-1,12),(-3,11),(-6,9),(-7,8),(-8,6),(-8,5),(-7,3),(-6,2),(-4,1),(-1,0)]
• 23-gon, area=387/2: [(0,0),(3,1),(5,2),(6,3),(7,5),(8,8),(8,9),(7,11),(6,12),(4,13),(1,14),(0,14),(-4,13),(-7,12),(-9,11),(-10,10),(-11,8),(-11,7),(-10,5),(-9,4),(-6,2),(-4,1),(-1,0)]
• 24-gon, area=420/2: [(0,0),(3,1),(5,2),(8,4),(9,5),(11,8),(12,10),(13,13),(13,14),(12,16),(11,17),(9,18),(8,18),(5,17),(3,16),(0,14),(-1,13),(-3,10),(-4,8),(-5,5),(-5,4),(-4,2),(-3,1),(-1,0)]
• 25-gon, area=497/2: [(0,0),(3,1),(5,2),(8,4),(9,5),(10,7),(11,10),(11,11),(10,13),(9,14),(7,15),(4,16),(3,16),(-1,15),(-4,14),(-6,13),(-9,11),(-10,10),(-11,8),(-11,7),(-10,5),(-9,4),(-6,2),(-4,1),(-1,0)]
• 26-gon, area=548/2: [(0,0),(3,1),(5,2),(8,4),(9,5),(11,8),(12,10),(13,13),(13,14),(12,16),(11,17),(9,18),(6,19),(5,19),(2,18),(0,17),(-3,15),(-4,14),(-6,11),(-7,9),(-8,6),(-8,5),(-7,3),(-6,2),(-4,1),(-1,0)]
• 27-gon, area=625/2: [(0,0),(3,1),(5,2),(8,4),(9,5),(11,8),(12,10),(13,13),(14,17),(14,18),(13,20),(12,21),(10,22),(9,22),(6,21),(4,20),(1,18),(-3,15),(-4,14),(-6,11),(-7,9),(-8,6),(-8,5),(-7,3),(-6,2),(-4,1),(-1,0)]
• 28-gon, area=690/2: [(0,0),(4,1),(7,2),(9,3),(12,5),(13,6),(15,9),(16,11),(17,14),(17,15),(16,17),(15,18),(13,19),(10,20),(9,20),(5,19),(2,18),(0,17),(-3,15),(-4,14),(-6,11),(-7,9),(-8,6),(-8,5),(-7,3),(-6,2),(-4,1),(-1,0)]
• 29-gon, area=783/2: [(0,0),(3,1),(5,2),(8,4),(9,5),(12,9),(14,12),(15,14),(16,17),(16,18),(15,20),(14,21),(12,22),(9,23),(8,23),(4,22),(1,21),(-1,20),(-4,18),(-5,17),(-7,14),(-8,12),(-9,9),(-9,8),(-8,5),(-7,3),(-6,2),(-4,1),(-1,0)]
• 30-gon, area=860/2: [(0,0),(4,1),(7,2),(9,3),(12,5),(13,6),(15,9),(16,11),(17,14),(17,15),(16,18),(15,20),(14,21),(12,22),(9,23),(8,23),(4,22),(1,21),(-1,20),(-4,18),(-5,17),(-7,14),(-8,12),(-9,9),(-9,8),(-8,5),(-7,3),(-6,2),(-4,1),(-1,0)]
• 31-gon, area=967/2: [(0,0),(4,1),(7,2),(9,3),(12,5),(13,6),(15,9),(16,11),(17,14),(18,18),(18,19),(17,22),(16,24),(15,25),(13,26),(12,26),(8,25),(5,24),(3,23),(0,21),(-4,18),(-5,17),(-7,14),(-8,12),(-9,9),(-9,8),(-8,5),(-7,3),(-6,2),(-4,1),(-1,0)]
• 32-gon, area=1046/2: [(0,0),(3,1),(5,2),(8,4),(9,5),(11,8),(12,10),(13,13),(13,14),(12,17),(11,19),(9,22),(8,23),(5,25),(3,26),(0,27),(-1,27),(-4,26),(-6,25),(-9,23),(-10,22),(-12,19),(-13,17),(-14,14),(-14,13),(-13,10),(-12,8),(-10,5),(-9,4),(-6,2),(-4,1),(-1,0)]
• 33-gon, area=1177/2: [(0,0),(3,1),(5,2),(8,4),(12,7),(13,8),(16,12),(18,15),(19,17),(20,20),(21,24),(21,25),(20,27),(19,28),(17,29),(14,30),(13,30),(9,29),(6,28),(4,27),(1,25),(-3,22),(-4,21),(-6,18),(-7,16),(-8,13),(-9,9),(-9,8),(-8,5),(-7,3),(-6,2),(-4,1),(-1,0)]
• 34-gon, area=1264/2: [(0,0),(4,1),(7,2),(9,3),(12,5),(13,6),(15,9),(16,11),(17,14),(17,15),(16,18),(15,20),(13,23),(12,24),(9,26),(7,27),(4,28),(3,28),(-1,27),(-4,26),(-6,25),(-9,23),(-10,22),(-12,19),(-13,17),(-14,14),(-14,13),(-13,10),(-12,8),(-10,5),(-9,4),(-6,2),(-4,1),(-1,0)]
• 35-gon, area=1409/2: [(0,0),(4,1),(7,2),(9,3),(12,5),(13,6),(15,9),(16,11),(17,14),(17,15),(16,18),(15,20),(13,23),(12,24),(10,25),(7,26),(3,27),(2,27),(-3,26),(-7,25),(-10,24),(-12,23),(-15,21),(-16,20),(-17,18),(-18,15),(-18,14),(-17,11),(-16,9),(-14,6),(-13,5),(-10,3),(-8,2),(-5,1),(-1,0)]
• 36-gon, area=1498/2: [(0,0),(4,1),(7,2),(9,3),(12,5),(13,6),(15,9),(16,11),(17,14),(17,15),(16,18),(15,20),(13,23),(12,24),(9,26),(7,27),(4,28),(0,29),(-1,29),(-5,28),(-8,27),(-10,26),(-13,24),(-14,23),(-16,20),(-17,18),(-18,15),(-18,14),(-17,11),(-16,9),(-14,6),(-13,5),(-10,3),(-8,2),(-5,1),(-1,0)]
• 37-gon, area=1671/2: [(0,0),(4,1),(7,2),(9,3),(12,5),(13,6),(16,10),(18,13),(19,15),(20,18),(20,19),(19,22),(18,24),(17,25),(14,27),(12,28),(9,29),(5,30),(4,30),(-1,29),(-5,28),(-8,27),(-10,26),(-13,24),(-14,23),(-16,20),(-17,18),(-18,15),(-18,14),(-17,11),(-16,9),(-14,6),(-13,5),(-10,3),(-8,2),(-5,1),(-1,0)]
• 38-gon, area=1780/2: [(0,0),(4,1),(7,2),(9,3),(12,5),(13,6),(15,9),(16,11),(17,14),(18,18),(18,19),(17,22),(16,24),(14,27),(13,28),(10,30),(8,31),(5,32),(1,33),(0,33),(-4,32),(-7,31),(-9,30),(-12,28),(-13,27),(-15,24),(-16,22),(-17,19),(-18,15),(-18,14),(-17,11),(-16,9),(-14,6),(-13,5),(-10,3),(-8,2),(-5,1),(-1,0)]
• 39-gon, area=1955/2: [(0,0),(4,1),(7,2),(9,3),(12,5),(13,6),(15,9),(16,11),(17,14),(18,18),(18,19),(17,22),(16,24),(14,27),(13,28),(10,30),(8,31),(5,32),(1,33),(0,33),(-5,32),(-9,31),(-12,30),(-14,29),(-17,27),(-18,26),(-20,23),(-21,21),(-22,18),(-22,17),(-21,14),(-20,12),(-18,9),(-17,8),(-13,5),(-10,3),(-8,2),(-5,1),(-1,0)]
• 40-gon, area=2078/2: [(0,0),(4,1),(7,2),(9,3),(12,5),(13,6),(15,9),(16,11),(17,14),(18,18),(18,19),(17,23),(16,26),(15,28),(13,31),(12,32),(9,34),(7,35),(4,36),(0,37),(-1,37),(-5,36),(-8,35),(-10,34),(-13,32),(-14,31),(-16,28),(-17,26),(-18,23),(-19,19),(-19,18),(-18,14),(-17,11),(-16,9),(-14,6),(-13,5),(-10,3),(-8,2),(-5,1),(-1,0)]

The target value of N can be given on the command line.

The minimum area of a convex polygon of given lattice width

The correctness of the program (knowing when one can stop) depends crucially on two lemmas that relate area to lattice width. The first lemma is about centrally symmetric polygons (arising when n is even), and the second lemma about arbitrary polygons.

Definition: The lattice width of a lattice polygon P in a given lattice direction (a,b) in Z² is one less than the number of lattice lines orthogonal to (a,b) intersecting P. The lattice width of P is the smallest number lattice width over all lattice directions.

Lemma 3:

A centrally symmetric convex polygon P of lattice width w has area at least w²/2.

Proof: [ A weaker bound of w²/4 is given in the following paper (p.175, before remark 2):

• Imre Bárány and Norihide Tokushige, The minimum area of convex lattice n-gons. Combinatorica 24 (No. 2, 2004), 171-185. https://www.renyi.hu/~barany/cikkek/94.pdf.

Our proof uses the same setup but refines the argument in the end. The improvement to w²/2 speeds up the computations by a factor 5-6. ]

1. We assume that the lattice width direction is vertical. Thus P touches two horizontal lines H+: y=b and H-: y=-b, with b=w/2. Let (a,0) be the intersection of P with the positive x-axis.

2. If a>=b, then the area is at least 2b²=w²/2, and we are done. Thus, let us assume a<b.

3. By horizontal shearing, we assume that the tangent at (a,0) (or some supporting line) is not vertical and has slope s at least +1. (Since horizontal lattice-preserving shearing changes the "inverse slope" 1/s in increments of 1, the inverse slope can always be brought into the interval (0,1].)

4. The lattice width in direction (1,0) is b'>=b. Thus P touches two vertical lines V+: x=b' and V-: x=-b'

5. The lattice width in direction (1,-1) is b">=b. Thus P touches two diagonal lines D+: y-x=b" and D-: y-x=-b" (Draw a picture!)

6. If (a,0) had tangent of slope 1, P could not touch D-, since a<b".

   Thus the slope s is strictly larger than 1, and P must touch V+ in the upper half-plane, and P must touch D- in the lower half-plane.

7. If the intersection of V+ and D- is below the x-axis, we swap the roles of V+ and D- by a horizontal shearing transformation combined with a vertical reflection: Each point (x,y) is mapped to (x-y,-y). This turns V+ into a line of slope 1 and makes D- vertical. Point on the x-axis are fixed.

   Thus we can assume that the intersection of V+ and D- is above or on the x-axis, at (b',c) with c>=0.

8. Pick a point where P touches each of the six lines, in a symmetric way. Denote them by B,BR,R (Bottom, Bottom-Right, and Right) and the reflected points by T,TL,L. (Draw a picture!) Replace P by the convex hull of these points (making the area smaller).

9. Holding L,B,R, and T fixed, we can move the points BR and TL on their respective lines D- and D+, keeping them symmetric. The area depends linearly on the movement.

   If the are would decrease when BR is moved towards R and towards V+, we can stop this movement at the point (b",0) where D- intersects the x-axis. The convex hull of (b",0),(-b",0),T,B has area 2b"b>w²/2.

10. Thus we are left with the case that we can decrease the area by moving BR towards H-, reaching the point B'=(d,-b) with d := b"-1>=0. Then we consider the convex hull of B', R, and their two symmetric points. R has coordinates R=(b',e) with e>=c>=0. The determinant of (d,-b) and (b',e) is bb'+de>=bb'>=b², and we are done.

QED.

Lemma 4:

A (not necessarily symmetric) convex polygon P of lattice width w has area at least w²/3. The inequality is strict except when P is a triangle.

Proof: The polygon P and the symmetrized version Q = 1/2 * (P-P) (which is half the "difference body" P-P={ a-b | a,b in P }), have the same width in every direction, and Q is centrally symmetric. We can apply Lemma 3 to Q. The area of Q is at most 3/2 times the area of P, and equality holds only when P is a triangle. See any of the following papers. QED

• Hans Rademacher, Über den Vektorenbereich eines konvexen ebenen Bereiches. Jahresbericht d. Deutschen Math.-Verein. 34 (1926): 64-78. http://eudml.org/doc/145705
• Theodor Estermann. Zwei neue Beweise eines Satzes von Blaschke und Rademacher. Jahresbericht d. Deutschen Math.-Verein. 36 (1927): 197-200. http://eudml.org/doc/145763
• Theodor Estermann. Über den Vektorenbereich eines konvexen Körpers. Math Z. 28 (1928), 471-475. https://doi.org/10.1007/BF01181177
• C. A. Rogers and G. C. Shephard, The difference body of a convex body. Arch. Math. 8 (1957), 220-233. https://doi.org/10.1007/BF01899997

I believe the true factor for Lemma 4 might be 3/8, and it is sharp for the triangle (-w/2,-w/2),(w/2,0),(0,w/2), for the case of even w. (By analogy with the ordinary width, not lattice width: "It is well-known that of all convex sets of a given width, the equilaterial triangle has the smallest area.")

In the proof, we would consider the 8-gon formed as the intersection of a horizontal slab of width w, a vertical slab of width w1>=1, and two slabs of directions 45° and -45° of vertical widths w3>=w and w4>=w, respectively. (Their Euclidean widths are >=w/sqrt(2).)

Some sides of the 8-gon may degenerate into a point, in case several lines are concurrent. The body must touch each of these edges (even if they are degenerated to a point).

Possible further improvements

Lemma 5.

A) For three successive vertices v,v',v" of P, the triangle vv'v" contains no interior point.

B) If n is even, then there cannot even be points on the diagonal vv".

Proof: A) Otherwise v' could be replaced by that point, reducing the area.

B) Bárány and Tokushige proved that the set E of edge vectors is convex in the following sense: Every vector e in the convex hull of E that is a primitive vector, must belong to E. Let e=v->v' and e'=v'->v". If vv"=e+e' would have gcd d>2 the vector (e+e')/d, which is in the convex hull of e,e',-e,-e', is missing from E, a contradiction. QED

This Lemma could be used to shortcut the propagation process. Property B is really strong. It seems that it holds also in the odd case, but I cannot prove it.

Property B means the following. Consider a pair (q,f). Then the outgoing edges q->p that need to be considered end on the next lattice line parallel to f. The points q form an arithmetic progression p0, p0+f, p0+2f, .... When considering a range of diameter D (roughly the maximum height), then these are about D/|f| successor points altogether (on all levels py).

Under Property A, we have to consider in addition all multiples of the vectors (q-f)->(p0+tf). I.e., p = (q-f) + s(p0+tf-(q-f)) = sp0 - (s-1)q + (st-s-1)f This would be D/|f| log(D|f|) ~ D log D/|f| points. Typically, the average length of the lists (the number of vectors f for a point q) is very small, (about 4 for height py=100), and |f| is usually not very short.

(By contrast, currently, every point p gets O(D²) inputs, one from every other point q (looking at the propagation from the incoming side).)

Some further possible improvements are mentioned in the comments in the program.

Now here comes finally the program:

N_target = 44

"""
Lemma 1: Every edge is a primitive vector.

Lemma 2: If n is even, then P can be assumed to be centrally symmetric.

Lemma 3: A centrally symmetric convex polygon P of lattice width w has area
at least w²/2.

Lemma 4: A (not necessarily symmetric) convex polygon P of lattice width w
has area at least w²/3.
The inequality is strict except when P is a triangle.

"""
# Since we use "volume", which is TWICE the area, the following constants
# are twice as big as the constants in the lemmas.

lattice_width_factor_even = 1
lattice_width_factor_odd  = 2/3 + 0.00001
# +0.00001 is because of the strict inequality (except for triangles, but
# those are found in the first iteration.)

def lower_bound(k,w):
    "lower bound on vol (=2*area) for k-gon of lattice width w "
    if k%2:
        lattice_width_factor = lattice_width_factor_odd
    else:
        lattice_width_factor = lattice_width_factor_even
    LB = ceil( w**2 * lattice_width_factor )
    if k%2 != LB%2:
        LB +=1 # vol must be odd for odd k and even for even k
    return LB

"""
===================================================
Also in the algorithm, we make the assumption
that the lattice width direction is vertical.
(The smallest number of parallel lattice lines intersecting P
are the horizontal lines.)

We assume that the lowest point (or one of the two lowest points)
is the origin O=(0,0).


=====================================================

For each pair of points p and q and every k, we remember the smallest
counterclockwise k-gon (0,0), .... , q, p

Restriction:
* The successive points move strictly upwards.
  (Thus, p is the top point, and there is a "long edge" from O to p.)

Polygons with one or two horizontal edges will be treated specially.

By horizontal lattice-preserving shearing transformations, it suffices to
store results only for the range 0 <= px < py
(The answer for any point p in the plane can then be figured out
by an appropriate shearing transformation.)

The area is roughly n^3, and by Lemma 4 the maximum height H=py can be bounded by
H=n^1.5 (omitting all constant factors). Thus the number of points p is n^3.
The points q can be bounded only very loosely in the horizontal
direction, by n^2, however, in practice, it turns out they go only up to
(1+alpha)n, where alpha grows very slowly (less than 4 for height 120.)

Thus, the number of points q that are considered is also roughly of the order n^3.
This gives an estimated runtime of O(n^3*n^3*n)=O(n^7) for all triples (p,q,k).
One could say O(n^6) for each new value of n.

===========================

In the end we put together a counterclockwise k-gon ending in p
with a counterclockwise k'-gon starting in O or in (-1,0) and
ending in p or in p-(1,0).

Main data structure:
====================
min_gon is a dictionary
min_gon[(px,py)] is a dictionary

min_gon[(px,py)][k] is a list of pairs (f,vol)=((fx,fy),vol), where
the incoming edge f=(fx,fy)=q-p is a primitive vector (pointing upwards, fy>0)
and vol = 2 * min-area of a k-gon ((0,0), .... ,q,p).
In the list, the vectors f are sorted clockwise by direction,
and the corresponding values "vol" are strictly increasing.
(If pairs don't fit this order, we can eliminate one of them, as DOMINATED.)

We fill this dictionary row by row, increasing py.

[ If we want to enumerate ALL optimal solutions, we should allow
the "vol" values to be weakly increasing. We should then use the
weaker lower bound also for even k, in case we are interested in solutions
that are not centrally symmetric.
Another possibility is COUNTING the number of solutions. (For each height
separately. Every solution has up to 2k representations with at least
one flat edge, and infinitely many without flat edge.)
We should consider all solutions of height up until the value is confirmed. ]
[ We should have proceeded from left to right, swapping x and y.
Then the notion of slope could have been used more naturally. ]

"vol" is TWICE the area.

https://oeis.org/A070911
1,2,5,6,13,14,21,28,43,48,65,80,103,118, / 151,174,213,242,289,328,387,420,497
The odd values after the slash were unconfirmed according to OEIS.
"""

known_value = {i+3:v for i,v in enumerate((
    1,2,5,6,13,14,21,28,43,48,65,80,103,118,151,174,213,242,289,
    328,387,420,497,548,625,690,783,860,967,1046,1177,1264,1409,
    ))} # just for guidance; not used by the algorithm

from math import tan,pi,gcd,sqrt,ceil
import itertools
import sys

if len(sys.argv)>1:
    N_target = int(sys.argv[1])

min_gon = dict()
record_vol = [n**3 for n in range(N_target+1)] # loose upper bound as start value
min_height = [0 for n in range(N_target+1)] # smallest height of a min-area polygon
max_height = [0 for n in range(N_target+1)]
# largest height of a min-area polygon with at least one horizontal edge,
# AS FOUND BY THE PROGRAM, not definite

confirmed = [True]*2 + [False]*(N_target+1-2) # checks which entries are confirmed.
how_achieved = [None]*(N_target+1)

# Auxiliary procedures to construct the solution, once the optimal
# area has been determined
def find_polygon(px,py,k,alpha0=0):
    # find the smallest k-gon ending in (px,py) by backtracing
    # 0<px<py. The two points (0,0),(1,0) are not produced.
    alphaT = alpha0
    result = [(px+alphaT*py,py)]
    (fx,fy),vol = min_gon[px,py][k][0]
    for k in range (k-1,1,-1):
        qx,qy = px-fx,py-fy
        vol -= qx*py-qy*px
        alpha =  qx//qy
        px,py = qx-alpha*qy,qy
        alphaT += alpha # alpha is accumulated
        result.append((px+alphaT*py,py))
        for (fx,fy),vol2 in min_gon[px,py][k]:
            if vol2==vol: # this must be the right entry, recognition by vol is easiest
                break
        else:
            error("not found",k,px,py)
    return result

def print_solution(k,how_achieved,vol):
    how,px,py,k1,k2 = how_achieved
    #print(f"{(k,how,px,py,k1,k2)=}")
    p1 = find_polygon(px,py,k1)
    p1.reverse()
    if how=="1 FLAT":
        top_point = px-1
    else:
        top_point = px
    ppx = (-top_point)%py
    alpha0 = (-top_point - ppx) // py # usually, alpha0==-1
    leftshift = 1 if how=="2 FLAT" else 0
    p2 = [(-x-leftshift,y) for x,y in find_polygon(ppx,py,k2,alpha0)]
    #print(f"-- {how=} {alpha0=}, {p2[0][0]=},{top_point=},{ppx=},{px=},{py=},{k=}")
    assert p2[0][0]==top_point-leftshift
    if how=="DIAGONAL":
        p2 = p2[1:] #eliminate the common point
    last = [(-1,0)] if how=="2 FLAT" else []
    print ("--- a smallest %d-gon,"%k, "area = %d/2:"%vol,
           [(0,0)]+p1+p2+last) # counterclockwise

## Start the computation

for py in itertools.count(1): # for py = 1,2,...
    # proceed row by row
    for px in range(py): # 2-gons are the base case.
        min_gon[px,py]={2:[((px,py),0)]}
        # or perhaps rather start 2-gons?

    work = num_work = 0
    sum_alpha = num_alpha = 0
    visited_p = visited_q = 0
    for k in range(2,N_target):
        collect=[[] for px in range(py)]
        # We will first *collect* the contributions to row py
        # for k+1 before sorting and processing them.         
        for qy in range(1,py):
            # Generate all contributions from point q on row qy to points p on row py.
            # We look at each starting point (qx0,qy) with 0<=qx0<qy<py, for which we
            # have min_gon[qx0,qy] stored.
            # Then we consider the series of affine images
            #    (qx0,qy), (qx0+qy,qy), (qx0+2*qy,qy), ..., (qx0+alpha*qy,qy), ...
            # until we are sure that there are no more contributions to row py
            # to be expected.
            #
            # These are two crucial pieces of the algorithm:
            # 1. The simultaneous clockwise sweep around each point (qx,qy)
            # with the incoming directions (fx,fy) and with the outgoing edges
            # to the points on row py.
            # 2. The logic to determine when no more progress can be expected from
            # advancing alpha.
            #            
            Delta_y = py-qy # >0
            for qx0 in range(qy):
                if k not in min_gon[qx0,qy]:
                    continue
                sides = min_gon[qx0,qy][k]
                assert(sides)
                for alpha in itertools.count(0): # for alpha = 0,1,2,...
                    qx = qx0+alpha*qy # apply shearing by alpha
                    visited_q += 1 # statistics
                    
                    px = 0
                    row_py_finished = False # for exit from nested loops
                    exhausted = True
                    for (fx0,fy),vol in sides:
                        # sweep clockwise around q; vol is increasing.
                        # simultaneously sweep the point p on
                        # row py from left to right, clockwise around q.
                        fx = fx0+alpha*fy # apply the reverse shearing by -alpha
                        while (# vector q->p is counterclockwise from (fx,fy):
                                # Delta_x/Delta_y < fx/fy 
                                (px-qx) * fy < fx * Delta_y ):
                            # store this:
                            Delta_x = px-qx
                            if gcd(Delta_x,Delta_y)==1:
                                collect[px].append(
                                    (vol + qx*py-qy*px, Delta_x,Delta_y))
                            if px>=py-1: # py-1 is largest value of px
                                row_py_finished = True
                                break
                            px += 1
                        if row_py_finished: break
                        
                        # Not all points on row py could be reached
                        # from the first incoming vector in q:
                        exhausted = False
                        
                    if exhausted:
                        # All points up to the maximum px have been reached
                        # by extending the FIRST incoming edge in "sides".
                        # Further increase of alpha would only produce a larger
                        # area of Opq, and could not lead to an improvement.
                        sum_alpha += alpha # statistics
                        num_alpha += 1
                        break
                    # [ There might be additional possibilities of shortcutting:
                    # For example, the volume of Opq is too big in
                    # order to be useful for anything. Indeed the lists in "collect"
                    # accumulate thousands of elements (from every point q),
                    # of which only 2-3 remain on average. ]

        # Now, consolidate the lists
        for px,triples in enumerate(collect):
            if not triples:
                continue
            visited_p += 1
            work += len(triples)
            num_work += 1 # statistics

            # This is one more crucial piece of the program:
            # We should keep only those incoming edges that are not DOMINATED.
            # If (fx',fy') enters clockwise from (fx,fy) and has larger volume,
            # then it cannot be part of a minimum-area solution.
            # If it has equal volume, it could be useful, but in any situation it
            # could always be replaced by (fx,fy).
            #
            # We solve this but looking at the entries sorted by slope.
            #
            result = []
            prev_fx, prev_fy = -1,0 # ensure that test is true at the first iteration
            prev_vol = -1
            for vol,fx,fy in (sorted(triples)):
                # smallest area first <==> result should be ordered clockwise.
                if (# this is the first time OR (fx,fy) is clockwise from prev:
                    # fx/fy > prev_fx/prev_fy
                      fx*prev_fy > prev_fx*fy):
                    if vol!=prev_vol:                   
                        result.append(((fx,fy),vol))
                    else: # replace the last entry:
                        result[-1] = ((fx,fy),vol)
                        # adjust this treatment when generating ALL solutions:
                        # Then multiple vol-entries should be kept, but SORTED BY SLOPE.
                    prev_fx,prev_fy,prev_vol = fx,fy,vol
                    
            min_gon[px,py][k+1]=result
            if 0:
                #print(f"p=({px},{py}) {k=}")
                for t in result: print("  ",t)

    # compute new record areas
    improved = [False]*(N_target+1)

    def check_record(k, vol, how, text=""):        
        if k>N_target:
            return
        if vol<record_vol[k]:
            #if text: print(text)
            assert not confirmed[k]

            if k in known_value and vol<known_value[k]:
                # then something must be wrong:
                print("BETTER THAN KNOWN VALUE",100*"!!")
                print_solution(k,how,vol)
                exit()              
            
            record_vol[k] = vol
            how_achieved[k] = how
            min_height[k] = py
            improved[k] = True
        if vol==record_vol[k] and how[0] != "DIAGONAL":
            # DIAGONAL case not counted: Without a flat edge, height is unbounded.
            max_height[k] = py

    # Combine "right" k1-gons with "left" k2-gons
    for px in range(py):
        ppx0 = (-px)%py  # A) exact match along long edge (0,0)-(px,py)
        ppx1 = (1-px)%py # B) horizontal edge of length 1 at the top
        for k1,best1 in min_gon[px,py].items():
            _,vol1 = best1[0]
            for k2,best2 in min_gon[ppx0,py].items():
                _,vol2 = best2[0]
                check_record(k1+k2-2,
                             vol1+vol2, ("DIAGONAL",px,py,k1,k2))
                check_record(k1+k2,
                             vol1+vol2+2*py, ("2 FLAT",px,py,k1,k2))
                 #f"DIAG k={k1+k2-3} {k1=} {k2=} {vol1=} {vol2=} p=({px},{py})")
            for k2,best2 in min_gon[ppx1,py].items():
                _,vol2 = best2[0]
                check_record(k1+k2-1, vol1+vol2+py, ("1 FLAT",px,py,k1,k2))
                 #f"EDGE k={k1+k2-2} {k1=} {k2=} {vol1=} {vol2=} p=({px},{py})")

    target_height = 0 # how high we still have to go
    for k in range(2,N_target+1):
        if not confirmed[k]:
            if lower_bound(k,py+1)>=record_vol[k]:
                # evaluated for increased py in the next iteration
                #print("LB",k,py+1,LB)
                confirmed[k] = True
                print_solution(k,how_achieved[k],record_vol[k])
            else:
                # find necessary target height to confirm the current record:
                t_k = py+2
                while lower_bound(k,t_k)<record_vol[k]:  t_k += 1
                target_height = max(target_height,t_k-1) 
                                
    print("N=%d, height %d"%(N_target,py) +
          (".\nRESULT:" if target_height==0 else "->%d:"%target_height),
          ",".join(("*" if improved[k] else "")+
                   str(record_vol[k])+
                   ("-" if k in known_value and record_vol[k]>known_value[k]
                    else "" if confirmed[k] else "?")
              for k in range(3,N_target+1)),
          #
          # Meaning of the signs:
          # - if there is a smaller known value
          # ? if still unconfirmed
          # * if has just been improved
          #
          "work=%d"%work, # work done in this step, total items in "collect"
          "" if num_work==0 else "avg=%.2f"%(work/num_work), # average list size in "collect"
          "" if num_alpha==0 else ("av_alpha=%.2f"% (sum_alpha/num_alpha)),
          flush=True)
    if py%10==0:
        print("smallest heights:",*("%2d"%x for x in min_height[3:]))
        print("largest heights: ",*("%2d"%x for x in max_height[3:]))
        print("2n-gons:", ",".join(str(record_vol[k]//2)
                                   +("" if confirmed[k] else "?")
                           for k in range(4,N_target+1,2)))

    if all(confirmed):
        break
    
print("smallest heights:",*("%2d"%x for x in min_height[3:]))
print("largest heights: ",*("%2d"%x for x in max_height[3:]))
# "largest heights" are preliminary, since the algorithm stops
# too early to find the largest heights.

# There is quite a large variation of heights, in line with the
# result of Bárány and Tukushige that, in the limit,
# the shapes resemble ellipses that are more and more oblong.

print("2n-gons:", ",".join(str(record_vol[k]//2)
                           for k in range(4,N_target+1,2)))

Answer 4

Rust + itertools + fasthash, n = 26 in ~423s

This is a rust port of Arnauld's answer, with a few iterations of performance improvement.

I use Vec<Option<_>> for sparse arrays and reference counting for the cache. I also use 128 * w + h as a hash for coordinate pairs, and use SeaHash from the fasthash crate

Obviously needs to be build with --release.

It is around 40 times faster for n=16, taking only 2.5s instead of 100s in Javascript.

use std::collections::{HashMap, HashSet};
use std::hash::{Hash, Hasher};
use std::iter::once;
use std::rc::Rc;
use std::time::{Duration, Instant};

use fasthash::sea::Hash64;
use itertools::Itertools;

struct Solution {
    w: i32,
    h: i32,
    y0: i32,
    res: [Vec<(i32, i32)>; 4],
}

#[derive(Clone, Debug)]
struct AreaPath {
    area: i32,
    path: Vec<(i32, i32)>,
}

#[derive(Eq, PartialEq)]
struct Coord(i32, i32);

#[allow(clippy::derive_hash_xor_eq)]
impl Hash for Coord {
    fn hash<H: Hasher>(&self, state: &mut H) {
        state.write_i32(128 * self.0 + self.1);
    }
}

type PathCache = HashMap<Coord, Rc<Vec<Option<AreaPath>>>, Hash64>;

fn main() {
    const MAX_N: i32 = 99;
    const MAX_TIME: Duration = Duration::from_secs(600);

    let mut pathCache: PathCache = HashMap::with_hasher(Hash64);

    let mut getPaths = |w: i32, h: i32| -> Rc<Vec<Option<AreaPath>>> {
        let cached: Option<&Rc<Vec<Option<AreaPath>>>> = pathCache.get(&Coord(w, h));
        return match cached {
            Some(p) => Rc::clone(p),
            None => {
                let mut set: HashSet<usize> = HashSet::new();
                let mut list = vec![None; 11];

                fn search(
                    x: i32,
                    y: i32,
                    pw: i32,
                    ph: i32,
                    area: i32,
                    path: Vec<(i32, i32)>,
                    w: i32,
                    h: i32,
                    set: &mut HashSet<usize>,
                    list: &mut Vec<Option<AreaPath>>,
                ) {
                    let len: usize = path.len();
                    if x == w && y == h {
                        let outerArea = w * h * 2 - area;

                        if !set.contains(&len) || outerArea > list[len].as_ref().unwrap().area {
                            set.insert(len);

                            list[len] = Some(AreaPath {
                                area: outerArea,
                                path: path.clone(),
                            });
                        }
                    }

                    for w0 in 1..=(w - x) {
                        for h0 in 1..=(h - y) {
                            if h0 * pw < w0 * ph {
                                let mut path_clone = Vec::with_capacity(len + 1);
                                path_clone.extend(path.clone());
                                path_clone.push((w0, h0));
                                search(
                                    x + w0,
                                    y + h0,
                                    w0,
                                    h0,
                                    area + ((w - x) * 2 - w0) * h0,
                                    path_clone,
                                    w,
                                    h,
                                    set,
                                    list,
                                );
                            }
                        }
                    }
                }

                search(0, 0, 0, 1, 0, vec![], w, h, &mut set, &mut list);
                let list = Rc::new(list);
                pathCache.insert(Coord(w, h), Rc::clone(&list));
                pathCache.insert(
                    Coord(h, w),
                    Rc::new(
                        list.iter()
                            .map(|o| {
                                o.as_ref().map(|a| AreaPath {
                                    area: a.area,
                                    path: a.path.iter().rev().map(|&(dx, dy)| (dy, dx)).collect(),
                                })
                            })
                            .collect(),
                    ),
                );

                list
            }
        };
    };

    let mut solveRectangle = |part: &[[i32; 4]],
                              w: i32,
                              h: i32,
                              y0: i32,
                              x0: i32,
                              y1: i32,
                              x1: i32|
     -> (f32, Option<[Vec<(i32, i32)>; 4]>) {
        let mut best = f32::INFINITY;
        let mut solution = None;

        for p in part {
            if let Some(Some(p0)) = getPaths(x0, h - y0 - 2).get((p[0] - 1) as usize) {
                if let Some(Some(p1)) = getPaths(w - x0 - 2, h - y1 - 1)
                    .get((p[1] - 1) as usize)
                {
                    if let Some(Some(p2)) = getPaths(w - x1 - 1, y1 - 1)
                        .get((p[2] - 1) as usize)
                    {
                        if let Some(Some(p3)) =
                            getPaths(x1 - 1, y0).get((p[3] - 1) as usize)
                        {
                            let outerArea = p0.area + p1.area + p2.area + p3.area;
                            let score = (w - 1) * (h - 1) * 2 - outerArea;
                            let score = score as f32;

                            if score < best {
                                best = score;
                                solution = Some([
                                    p0.path.clone(),
                                    p1.path.iter().rev().copied().collect(),
                                    p2.path.clone(),
                                    p3.path.iter().rev().copied().collect(),
                                ]);
                            }
                        }
                    }
                }
            }
        }
        return (best, solution);
    };

    let mut solve = |n: i32, min: i32, max: i32| -> (f32, Solution) {
        let part: Vec<[i32; 4]> = partitionsInFour(n);
        let mut rect = vec![];
        let mut best = f32::INFINITY;
        let mut solution = None;

        for w in min..=max {
            for h in min..=w {
                rect.push((w, h));
            }
        }

        for (w, h) in rect {
            for y0 in 0..(h - 1) {
                for x0 in 0..(w - 1) {
                    for y1 in 1..h {
                        for x1 in 1..w {
                            let (score, res) = solveRectangle(&part, w, h, y0, x0, y1, x1);
                            if score < best {
                                best = score;
                                solution = Some(Solution {
                                    w: w,
                                    h: h,
                                    y0: y0,
                                    res: res.unwrap(),
                                });
                            }
                        }
                    }
                }
            }
        }
        return (best, solution.unwrap());
    };

    let ts = Instant::now();
    let mut min = 4;
    let mut max = 4;

    for n in 7..=MAX_N {
        println!("n = {n}");

        let (score, solution) = solve(n, min, max);
        let w = solution.w - 1;
        let h = solution.h - 1;

        println!("Area * 2 = {score}");
        println!("Enclosing rectangle: {w} x {h}");
        displaySolution(&solution);

        let time = ts.elapsed();
        println!("Total time: {time:?}");

        if time >= MAX_TIME {
            break;
        }
        min = solution.h;
        max = solution.w + 2;
    }

    fn partitionsInFour(n: i32) -> Vec<[i32; 4]> {
        let mut list = vec![];
        fn search(n: i32, i: i32, l: Vec<i32>, list: &mut Vec<[i32; 4]>) {
            if n != 0 {
                if i <= n && l.len() != 4 {
                    search(n - i, 1, once(i).chain(l.iter().copied()).collect(), list);
                    search(n, i + 1, l, list);
                }
            } else if let Ok(a) = l.try_into() {
                list.push(a);
            }
        }

        search(n, 1, vec![], &mut list);
        return list;
    }

    fn displaySolution(solution: &Solution) {
        let mut m = vec![vec!["--".to_owned(); solution.w as usize]; solution.h as usize];
        let mut x = 0;
        let mut y;
        let mut i = 1;

        let mut mark = |x: i32, y: i32| {
            let x = x as usize;
            let y = y as usize;
            if m[y][x] == "--" {
                m[y][x] = format!("{:>02}", i);
            }
            i += 1;
        };

        y = solution.y0;
        mark(0, y);
        y += 1;
        mark(x, y);
        for (dx, dy) in &solution.res[0] {
            x += dx;
            y += dy;
            mark(x, y)
        }
        x += 1;
        mark(x, y);
        for (dx, dy) in &solution.res[1] {
            x += dx;
            y -= dy;
            mark(x, y)
        }
        y -= 1;
        mark(x, y);
        for (dx, dy) in &solution.res[2] {
            x -= dx;
            y -= dy;
            mark(x, y)
        }
        x -= 1;
        mark(x, y);
        for (dx, dy) in &solution.res[3] {
            x -= dx;
            y += dy;
            mark(x, y)
        }

        println!("{}", m.into_iter().map(|r| r.join(" ")).join("\n"));
    }
}

Answer 5

Python (with restricted conditions) n=26, 10 seconds

Uses backtracking.
Numba for speed.
notebook link

10.49 seconds, 548 area
points: [[0, 0], [0, 1], [1, 3], [2, 4], [4, 5], [7, 6], [8, 6], [11, 5], [13, 4], [16, 2], [17, 1], [19, -2], [20, -4], [21, -7], [21, -8], [20, -10], [19, -11], [17, -12], [14, -13], [13, -13], [10, -12], [8, -11], [5, -9], [4, -8], [2, -5], [1, -3]]

Fast O(n^2) upper bounds for high n (<10% error for n<=36, slightly better than n^3/27 when n<1e5, sometimes better for n>1e5)

Assuming the solution to be close to circular. We can make circle-like shapes where sides made of the smallest n vectors where i and j are coprime. This gives a polygon of size n*4. I think this should have the optimal minimum perimeter as well. github notebook link

The calculated upper bound is (slightly) better than n^3/27 for n upto 10000.

(Old extra method - best solution found in 15s for n upto 44 link . At these high numbers, it's basically just oneshot generations but still somewhat interesting)

Old code links n = 14, 80s

---

88.12 seconds, 80 area, points: [[0, 0], [0, 1], [1, 2], [3, 3], [4, 3], [6, 2], [7, 1], [8, -1], [8, -2], [7, -3], [5, -4], [4, -4], [2, -3], [1, -2]]

notebook link

---

another way to get it to generate 1 valid, pretty optimal output very early is to sort the points in order of distance to previous point, but it doesn't help in overall speed. link

---

if we assume point n+1 must be in a 7x7 grid centered at point n, the max manhattan distance between points is <6 and that the angle between formed by n-1, n and n+1 is >= 90, we get the output for n=18, gridsize=14 in 360 seconds

362.6 seconds, 174 area, points: [[0, 0], [0, 1], [1, 3], [2, 4], [4, 5], [5, 5], [7, 4], [8, 3], [9, 1], [10, -2], [10, -3], [9, -5], [8, -6], [6, -7], [5, -7], [3, -6], [2, -5], [1, -3]]

notebook link

---

reducing the number of wasted calculations and function calls since cross product is the area of a triangle and taking into account that p0 = 0,0 we get:

from numba import jit
import numpy as np

@jit
def create_polygon_jit(points, n, target_n, gridsize, curr_area, min_area, soln):
    if n >= 3:
        if points[1,0]==0 and points[n-1,0]==0:
            return False
        if (points[n-3][0]-points[n-2][0])*(points[n-1][1]-points[n-2][1])<=\
        (points[n-3][1]-points[n-2][1])*(points[n-1][0]-points[n-2][0]):
            return False
        add_area = points[n-1][0]*points[n-2][1]-points[n-1][1]*points[n-2][0]
        if add_area <= 0:
            return False
        curr_area += add_area
        if curr_area >= min_area[0]:
            return False
    if n == target_n:
        min_area[0] = curr_area
        soln[:] = points.copy()
        return False
    for i in range(gridsize):
        for j in range(-gridsize//2, gridsize//2):
            points[n] = (i,j)
            result = create_polygon_jit(points, n+1, target_n, gridsize, curr_area, min_area, soln)
    points[n] = (0,0)
    return False

which only does 4 seconds better at 84.4 sec notebook link

---

pypy 3 tio link

---

restricted conditions with quadrant/octant optimization suggested by @HugoPfoertner, n=19 works, and n=18 is done in 180 seconds

438.59 seconds, 213 area, points: [[0, 0], [0, 1], [1, 3], [2, 4], [4, 5], [5, 5], [8, 4], [10, 3], [11, 2], [12, 0], [13, -3], [13, -4], [12, -6], [11, -7], [10, -7], [7, -6], [5, -5], [2, -3], [1, -2]]

notebook link

Answer 6

C++ (gcc), n=16 ~2 seconds

C++ port of @arrmansa's optimized version.

#include <array>
#include <chrono>
#include <iostream>

struct point {
    int x;
    int y;
};

template<std::size_t N>
bool create_polygon_jit(std::array<point, N>& points, int n, int target_n, int gridsize, int curr_area, int& min_area, std::array<point, N>& soln)
{
    if (n >= 3)
    {
        int x1 = points[n - 3].x - points[n - 2].x, y1 = points[n - 3].y - points[n - 2].y;
        int x2 = points[n - 1].x - points[n - 2].x, y2 = points[n - 1].y - points[n - 2].y;
        if (x1 * y2 - x2 * y1 <= 0) // Angle > 180
            return false;
        if (x1 * x2 + y1 * y2 > 0)  // Angle >= 90
            return false;
        if (points[1].x == 0 and points[n - 1].x == 0)
            return true;

        int add_area = points[n - 1].x * points[n - 2].y - points[n - 1].y * points[n - 2].x;
        if (add_area <= 0)
            return true;

        curr_area += add_area;
        if (curr_area >= min_area)
            return true;

        if (n == target_n)
        {
            min_area = curr_area;
            soln = points;
            return true;
        }

        int min_i = std::max(0, points[n - 1].x - 3);
        int max_i = std::min(gridsize, points[n - 1].x + 4);

        int min_j = std::max(-gridsize / 2, points[n - 1].y - 3);
        int max_j = std::min(gridsize / 2 + 1, points[n - 1].y + 4);

        if (x2 > 0)
        {
            if (y2 > 0)
            {
                if (x2 > y2)
                {
                    for (int j = min_j; j < max_j; j++) {
                        for (int i = max_i - 1; i > min_i - 1; i--) {
                            points[n].x = i;
                            points[n].y = j;
                            if (!create_polygon_jit(points, n + 1, target_n, gridsize, curr_area, min_area, soln))
                                break;
                        }
                        if (points[n].x == max_i - 1)
                            break;
                    }
                }
                else
                {
                    for (int j = max_j - 1; j > min_j - 1; j--)
                    {
                        for (int i = max_i - 1; i > min_i - 1; i--)
                        {
                            points[n].x = i;
                            points[n].y = j;

                            if (!create_polygon_jit(points, n + 1, target_n, gridsize, curr_area, min_area, soln))
                                break;
                        }
                        if (points[n].x == max_i - 1)
                            break;
                    }
                }
            }
            else
            {
                if (x2 > -y2)
                {
                    for (int i = max_i - 1; i > min_i - 1; i--)
                    {
                        for (int j = min_j; j < max_j; j++)
                        {
                            points[n].x = i;
                            points[n].y = j;

                            if (!create_polygon_jit(points, n + 1, target_n, gridsize, curr_area, min_area, soln))
                                break;
                        }
                        if (points[n].y == min_j)
                            break;
                    }
                }
                else
                {
                    for (int i = min_i; i < max_i; i++)
                    {
                        for (int j = min_j; j < max_j; j++)
                        {
                            points[n].x = i;
                            points[n].y = j;
                            ;
                            if (!create_polygon_jit(points, n + 1, target_n, gridsize, curr_area, min_area, soln))
                                break;
                        }
                        if (points[n].y == min_j)
                            break;
                    }
                }
            }
        }
        else
        {
            if (y2 > 0)
            {
                if (-x2 > y2)
                {
                    for (int i = min_i; i < max_i; i++)
                    {
                        for (int j = max_j - 1; j > min_j - 1; j--)
                        {
                            points[n].x = i;
                            points[n].y = j;

                            if (!create_polygon_jit(points, n + 1, target_n, gridsize, curr_area, min_area, soln))
                                break;
                        }
                        if (points[n].y == max_j - 1)
                            break;
                    }
                }
                else
                {
                    for (int i = max_i - 1; i > min_i - 1; i--)
                    {
                        for (int j = max_j - 1; j > min_j - 1; j--)
                        {
                            points[n].x = i;
                            points[n].y = j;

                            if (!create_polygon_jit(points, n + 1, target_n, gridsize, curr_area, min_area, soln))
                                break;
                        }
                        if (points[n].y == max_j - 1)
                            break;
                    }
                }
            }
            else
            {
                if (-x2 > -y2)
                {
                    for (int j = max_j - 1; j > min_j - 1; j--)
                    {
                        for (int i = min_i; i < max_i; i++)
                        {
                            points[n].x = i;
                            points[n].y = j;

                            if (!create_polygon_jit(points, n + 1, target_n, gridsize, curr_area, min_area, soln))
                                break;
                        }
                        if (points[n].x == min_i)
                            break;
                    }
                }
                else
                {
                    for (int j = min_j; j < max_j; j++)
                    {
                        for (int i = min_i; i < max_i; i++)
                        {
                            points[n].x = i;
                            points[n].y = j;

                            if (!create_polygon_jit(points, n + 1, target_n, gridsize, curr_area, min_area, soln))
                                break;
                        }
                        if (points[n].x == min_i)
                            break;
                    }
                }

            }
        }
    }
    else
    {
        int min_i = std::max(0, points[n - 1].x - 3);
        int max_i = std::min(gridsize, points[n - 1].x + 4);

        int min_j = std::max(-gridsize / 2, points[n - 1].y - 3);
        int max_j = std::min(gridsize / 2, points[n - 1].y + 4);

        for (int i = min_i; i < max_i; i++)
        {
            for (int j = min_j; j < max_j; j++)
            {
                points[n].x = i;
                points[n].y = j;
                create_polygon_jit(points, n + 1, target_n, gridsize, curr_area, min_area, soln);
            }
        }
    }
    points[n].x = 0;
    points[n].y = 0;
    return true;
}

int main() {
    using namespace std::chrono;
    auto t1 = steady_clock::now();

    const int n = 16;
    std::array<point, n> points = { 0 };

    int gridSize = 11;
    int minArea = gridSize * gridSize;

    std::array<point, n> result;


    create_polygon_jit(points, 1, n, gridSize, 0, minArea, result);

    std::cout << "Result: ";
    for (point& p : result)
        std::cout << "(" << p.x << ", " << p.y << ") ";
    std::cout << "\nArea: " << minArea << "\n";

    auto t2 = steady_clock::now();
    auto time_span = duration_cast<duration<double>>(t2 - t1);
    std::cout << "Time:" << time_span.count() << " seconds.\n";
}

Try it online!

Answer 7

Rust, n=44 in ~7 seconds

This is merely an abridged port of Günter Rote's solution, no original contribution.

const LWF_EVEN: f64 = 1.0;
const LWF_ODD: f64  = 2.0/3.0 + 0.00001;

fn lower_bound(k: usize, w: usize) -> usize {
    let lattice_width_factor = if k%2 == 1 {
        LWF_ODD
    } else {
        LWF_EVEN
    };
    let mut lb = ( (w*w) as f64 * lattice_width_factor ) as usize;
    if k%2 != lb%2 {
        lb += 1;
    };
    lb
}

fn gcd(mut x: isize, mut y: isize) -> isize {
    while y != 0 {
        let t = y;
        y = x % y;
        x = t;
    }
    x
}

fn check_record(k: usize, vol: usize, py: usize, n_target: usize,
    record_vol: &mut [usize], min_height: &mut [usize], max_height: &mut [usize], improved: &mut [bool], diagonal: bool) {
    if k > n_target {
        return
    }
    if vol < record_vol[k] {
        record_vol[k] = vol;
        min_height[k] = py;
        improved[k] = true;
    }
    if vol==record_vol[k] && !diagonal {
        max_height[k] = py;
    }
}

use std::collections::HashMap;

fn main() {
    #[allow(non_snake_case)]
    let mut N_target = 44;
    if let Some(s) = std::env::args().skip(1).next() {
        if let Ok(x) = s.parse() {
            N_target = x;
        }
    };
    let mut min_gon = HashMap::new();
    let mut record_vol = Vec::new();
    for n in 0..=N_target {
        record_vol.push(n*n*n);
    };
    let mut min_height = vec![0; N_target+1];
    let mut max_height = vec![0; N_target+1];
    let mut confirmed = vec![false; N_target+1];
    confirmed[0] = true; confirmed[1] = true;
    for py in 1usize.. {
        for px in 0usize..py {
            min_gon.insert((px, py), HashMap::from([(2, vec![((px as isize, py as isize), 0)])]));
        };

        let mut _work = 0;
        let mut _num_work = 0;
        let mut _sum_alpha = 0;
        let mut _num_alpha = 0;
        let mut _visited_p = 0;
        let mut _visited_q = 0;
        for k in 2usize..N_target {
            let mut collect = Vec::new();
            for _px in 0..py {
                collect.push(Vec::new());
            };
            for qy in 1..py {
                let delta_y = (py-qy) as isize;
                for qx0 in 0..qy {
                    let sides = match min_gon[&(qx0,qy)].get(&k) {
                        Some(x) => x,
                        None => continue
                    };
                    for alpha in 0.. {
                        let qx = qx0+alpha*qy;
                        _visited_q += 1;
                        let mut px = 0usize;
                        let mut exhausted = true;
                        'a: for ((fx0,fy),vol) in sides {
                            let fx = fx0+alpha as isize*fy;
                            while (px as isize - qx as isize) * fy < fx * delta_y {
                                let delta_x = px as isize - qx as isize;
                                if gcd(delta_x.abs(), delta_y) == 1 {
                                    collect[px].push(
                                        (vol + qx*py-qy*px, delta_x, delta_y));
                                };
                                if px>=py-1 {
                                    break 'a
                                }
                                px += 1;
                            }
                            exhausted = false;
                        }
                        if exhausted {
                            _sum_alpha += alpha;
                            _num_alpha += 1;
                            break
                        }
                    }
                }
            }
            for (px, mut triples) in collect.into_iter().enumerate() {
                if triples.is_empty() {
                    continue
                };
                _visited_p += 1;
                _work += triples.len();
                _num_work += 1;
                let mut result = Vec::new();
                let mut first_time = true;
                let (mut prev_fx, mut prev_fy) = (0, 0);
                let mut prev_vol = 0;
                triples.sort();
                for (vol,fx,fy) in triples.into_iter() {
                    if first_time || fx*prev_fy > prev_fx*fy {
                        if first_time || vol != prev_vol {
                            result.push(((fx,fy),vol));
                            first_time = false;
                        } else {
                            let i = result.len()-1;
                            result[i] = ((fx,fy),vol);
                        };
                        prev_fx = fx; prev_fy = fy; prev_vol = vol;
                    }
                };
                min_gon.get_mut(&(px,py)).expect("min_gon[px,py] has to be initialized").insert(k+1, result);
            }
        }
        let mut improved = vec![false; N_target+1];
        for px in 0..py {
            let ppx0 = (py - px) % py;
            let ppx1 = (py + 1 - px) % py;
            for (k1,best1) in min_gon[&(px,py)].iter() {
                let vol1 = best1[0].1;
                for (k2,best2) in min_gon[&(ppx0,py)].iter() {
                    let vol2 = best2[0].1;
                    check_record(k1+k2-2, vol1+vol2, py, N_target, &mut record_vol, &mut min_height, &mut max_height, &mut improved, true);
                    check_record(k1+k2, vol1+vol2+2*py, py, N_target, &mut record_vol, &mut min_height, &mut max_height, &mut improved, false);
                }
                for (k2,best2) in min_gon[&(ppx1,py)].iter() {
                    let vol2 = best2[0].1;
                    check_record(k1+k2-1, vol1+vol2+py, py, N_target, &mut record_vol, &mut min_height, &mut max_height, &mut improved, false);
                }
            }
        }
        let mut target_height = 0;
         for k in 2..=N_target {
            if !confirmed[k] {
                if lower_bound(k,py+1) >= record_vol[k] {
                    confirmed[k] = true;
                } else {
                    let mut t_k = py+2;
                    while lower_bound(k,t_k) < record_vol[k] { t_k += 1 };
                    target_height = std::cmp::max(target_height, t_k-1);
                }
            }
        }
        if confirmed.iter().all(|x| *x) {
            break
        }
    };
    println!("{:?}", &record_vol[3..]);
}