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..