Problem: https://open.kattis.com/problems/randommanhattan

Official Solution:

/*

Within a zone:  Z = (l^2+3lr+r^2)w^3/15

Cx(i) = x0 + (l+2r)/(3(l+r))w

*/
#include <iostream>
#include <algorithm>
#include <iomanip>
#include <cstdlib>
#include <vector>
#include <math.h>
#include <tuple>
using namespace std ;
using ll = long long ;
using t3 = tuple<double, double, double> ;
int main() {
   cout << setprecision(15) ;
   int N{0} ;
   cin >> N ;
   vector<ll> xs(N), ys(N) ;
   for (int i=0; i<N; i++)
      cin >> xs[i] >> ys[i] ;
   double r {} ;
   for (int outer=0; outer<2; outer++) {
      int hix0 = min_element(xs.begin(), xs.end()) - xs.begin() ;
      int lox0 { hix0 } ;
      auto ht=[&](ll x, int p0, int p1) -> double {
         if (x == xs[p0])
            return ys[p0] ;
         if (x == xs[p1])
            return ys[p1] ;
         return ys[p0]+(x-xs[p0])*(ys[p1]-ys[p0])/(double)(xs[p1]-xs[p0]) ;
      } ;
      vector<t3> zones ;
      while (1) {
         int hix1 { (hix0 + 1) % N } ;
         int lox1 { (lox0 + N - 1) % N } ;
         ll x0 = max(xs[hix0], xs[lox0]) ;
         if (x0 == xs[hix1]) {
            hix0 = hix1 ;
            continue ;
         }
         if (x0 == xs[lox1]) {
            lox0 = lox1 ;
            continue ;
         }
         ll x1 { min(xs[hix1], xs[lox1]) } ;
         if (x1 < x0)
            break ;
         if (x1 == x0)
            throw "Failed while building zones" ;
         double lft { ht(x0, hix0, hix1) - ht(x0, lox0, lox1) } ;
         double rgt { ht(x1, hix0, hix1) - ht(x1, lox0, lox1) } ;
         zones.push_back(make_tuple(lft, rgt, (double)(x1-x0))) ;
//  cout << "Adding tuple xs " << x0 << " " << x1 << " heights " << lft << " " << rgt << endl ;
         if (x1 == xs[hix1])
            hix0 = hix1 ;
         if (x1 == xs[lox1])
            lox0 = lox1 ;
      }
      double s {} ;
      double cxa {} ;
      double sa {} ;
      double a2 {} ;
      double x0 {} ;
      for (int i=0; i<zones.size(); i++) {
         t3 &z = zones[i] ;
         double lft { get<0>(z) }, rgt { get<1>(z) }, w { get<2>(z) } ;
         s += (lft*lft+3*lft*rgt+rgt*rgt)*w*w*w/15 ;
         double ta { (lft+rgt)*w/2 } ;
         double cx { x0 + (lft+2*rgt)/(3*(lft+rgt))*w } ;
         if (i)
            s += 2 * (cx - cxa / sa) * sa * ta ;
         a2 += ta*(2*sa+ta) ;
         sa += ta ;
         cxa += cx * ta ;
         x0 += w ;
      }
      r += s / a2 ;
      swap(xs, ys) ;
      reverse(xs.begin(), xs.end()) ;
      reverse(ys.begin(), ys.end()) ;
   }
   cout << r << endl ;
}

I understand that Cx is the center of a trapezoid but I don't know how the rest of the code works. I know they do 2 passes and swap x and y coordinates, but I don't know how they got the equation for Z (why the 15? does it involve integration?), why they're using the (current area) * (2 * previously accumulated area + current area), why they multiple the trapezoid center by the trapezoid area to get cxa, etc..