Solution -「多校联训」最大面积

\(\mathcal{Description}\)

  Link.

  平面上有 \(n\) 个点 \(A_{1..n}\)，\(q\) 次询问，每次给出点 \(P\)，求

\[\max_{1\le l\le r\le n}\left\{\sum_{i=l}^r \vec{OP}\times\vec{OA_i}\right\}. \]

  \(n\le10^5\)，\(q\le10^6\)。

\(\mathcal{Solution}\)

  初步转化一下式子：

\[\begin{aligned} \sum_{i=l}^r\vec{OP}\times\vec{OA_i}&=\sum_{i=l}^r(x_Py_{A_i}-y_Px_{A_i})\\ &=x_P\left(\sum_{i=l}^ry_{A_i}-\frac{y_P}{x_P}\sum_{i=l}^rx_{A_i}\right)~~~~(x_P\not=0) \end{aligned} \]

对于 \(x_P=0\) 即求 \(x_{A_{1..n}}\) 的最大子段和；对于 \(x_P\not=0\)，提出 \(x_P\) 后，记 \(k=\frac{y_P}{x_P}\)，\(y(l,r)=\sum_{i=l}^ry_{A_i}\)，\(x(l,r)=\sum_{i=l}^rx_{A_i}\)，原式可以看作

\[y(l,r)-kx(l,r)=b, \]

并要求最大化 \(x_Pb\)。相当于过某个 \((x(l,r),y(l,r))\) 作斜率为 \(k\) 的直线，求其最大或最小纵截距。若能求出所有 \((x(l,r),y(l,r))\) 构成的凸包，就能在凸壳上二分求解了。

  考虑求凸包的方法：分治，对于分支区间 \([l,r]\) 与其中点 \(p\)，计算 \((x(l..p,p+1..r),y(l..p,p+1..r))\) 的贡献。分别求出后缀与前缀凸包，根据定义，求出左右凸包的 Minkowski 和即为当前层贡献，加入总点集，最后再求总点集的凸包即可。

  复杂度 \(\mathcal O(n\log^2n+q\log n)\)。

\(\mathcal{Code}\)

/* Rainybunny */

#include <cstdio>
#include <vector>
#include <algorithm>

#define rep( i, l, r ) for ( int i = l, rep##i = r; i <= rep##i; ++i )
#define per( i, r, l ) for ( int i = r, per##i = l; i >= per##i; --i )

typedef long long LL;
typedef long double LD;

inline int rint() {
    int x = 0, f = 1, s = getchar();
    for ( ; s < '0' || '9' < s; s = getchar() ) f = s == '-' ? -f : f;
    for ( ; '0' <= s && s <= '9'; s = getchar() ) x = x * 10 + ( s ^ '0' );
    return x * f;
}

template<typename Tp>
inline void wint( Tp x ) {
    if ( x < 0 ) putchar( '-' ), x = -x;
    if ( 9 < x ) wint( x / 10 );
    putchar( x % 10 ^ '0' );
}

inline LL lmax( const LL a, const LL b ) { return a < b ? b : a; }

namespace PGP {

const LD EPS = 1e-9;

inline int dcmp( const LD a ) {
    return -EPS < a && a < EPS ? 0 : a < 0 ? -1 : 1;
}

struct Point {
    LL x, y;
    Point(): x( 0 ), y( 0 ) {}
    Point( const LL a, const LL b ): x( a ), y( b ) {}
    inline Point operator + ( const Point& p ) const {
        return { x + p.x, y + p.y };
    }
    inline Point operator - ( const Point& p ) const {
        return { x - p.x, y - p.y };
    }
    inline int operator ^ ( const Point& p ) const {
        return dcmp( LD( x ) * p.y - LD( y ) * p.x );
    }
    inline bool operator < ( const Point& p ) const {
        return x != p.x ? x < p.x : y < p.y;
    }
    inline bool operator == ( const Point& p ) const {
        return x == p.x && y == p.y;
    }
};
typedef Point Vector;
typedef std::vector<Point> Convex;

inline Convex getConvex( std::vector<Point> vec ) {
    static Convex ret; ret.clear();
    std::sort( vec.begin(), vec.end() );
    vec.resize( std::unique( vec.begin(), vec.end() ) - vec.begin() );
    int n = int( vec.size() ), top = 0;
    ret.resize( n << 1 );
    rep ( i, 0, n - 1 ) {
        for ( ; top > 1; --top ) {
            if ( ( ( ret[top - 1] - ret[top - 2] )
              ^ ( vec[i] - ret[top - 2] ) ) > 0 ) break;
        }
        ret[top++] = vec[i];
    }
    for ( int tmp = top, i = n - 2; ~i; --i ) {
        for ( ; top > tmp; --top ) {
            if ( ( ( ret[top - 1] - ret[top - 2] )
              ^ ( vec[i] - ret[top - 2] ) ) > 0 ) break;
        }
        ret[top++] = vec[i];
    }
    top -= n > 1;
    return ret.resize( top ), ret;
}

inline int findPole( const Convex& cvx ) {
    int ret = -1, n = int( cvx.size() );
    rep ( i, 0, n - 1 ) {
        if ( !~ret || cvx[ret].y > cvx[i].y
          || ( cvx[ret].y == cvx[i].y && cvx[ret].x > cvx[i].x ) ) {
            ret = i;
        }
    }
    return ret;
}

inline void poleSort( Convex& cvx ) {
    int pid = findPole( cvx ), n = int( cvx.size() );
    pid = n - pid - 1;
    std::reverse( cvx.begin(), cvx.end() );
    std::reverse( cvx.begin(), cvx.begin() + pid + 1 );
    std::reverse( cvx.begin() + pid + 1, cvx.end() );
}

inline Convex minkowskiSum( Convex A, Convex B ) {
    static Convex ret; ret.clear();
    poleSort( A ), poleSort( B );
    int n = int( A.size() ), m = int( B.size() );
    
    Point ap( A[0] ), bp( B[0] );
    ret.push_back( ap + bp );
    rep ( i, 0, n - 2 ) A[i] = A[i + 1] - A[i];
    A[n - 1] = ap - A[n - 1];
    rep ( i, 0, m - 2 ) B[i] = B[i + 1] - B[i];
    B[m - 1] = bp - B[m - 1];
    
    int i = 0, j = 0;
    while ( i < n && j < m ) {
        // 注意这里能 hack，应该计算极角来比较，否则无法正确处理共线向量。
        ret.push_back( ret.back()
          + ( ( A[i] ^ B[j] ) > 0 ? A[i++] : B[j++] ) );
    }
    while ( i < n ) ret.push_back( ret.back() + A[i++] );
    while ( j < m ) ret.push_back( ret.back() + B[j++] );
    return ret;
}

} using namespace PGP;

const int MAXN = 1e5;
int n, q, rpos;
Point A[MAXN + 5];
Convex cvxA;
LL mnsec, mxsec;

inline void buildConvex( const int l, const int r ) {
    if ( l == r ) return cvxA.push_back( A[l] );
    int mid = l + r >> 1;
    buildConvex( l, mid ), buildConvex( mid + 1, r );

    static Convex cvxL, cvxR; cvxL.clear(), cvxR.clear();
    Point cur;
    per ( i, mid, l ) cvxL.push_back( cur = cur + A[i] );
    cur = Point();
    rep ( i, mid + 1, r ) cvxR.push_back( cur = cur + A[i] );
    
    cvxL = minkowskiSum( getConvex( cvxL ), getConvex( cvxR ) );
    for ( const auto& p: cvxL ) cvxA.push_back( p );
}

inline void init() {
    buildConvex( 1, n ), cvxA = getConvex( cvxA );

    #ifdef RYBY
        puts( "+++ cvxA +++" );
        for ( auto p: cvxA ) printf( "%lld %lld\n", p.x, p.y );
        puts( "--- cvxA ---" );
    #endif

    if ( cvxA.size() == 1 ) rpos = 0;
    else {
        rpos = -1;
        rep ( i, 1, cvxA.size() - 1 ) {
            if ( cvxA[i].x <= cvxA[i - 1].x ) {
                rpos = i - 1; break;
            }
        }
        if ( !~rpos ) rpos = cvxA.size() - 1;
    }

    LL las = 0;
    rep ( i, 1, n ) {
        if ( ( las += A[i].x ) > mxsec ) mxsec = las;
        if ( las < 0 ) las = 0;
    }
    las = 0;
    rep ( i, 1, n ) {
        if ( ( las += A[i].x ) < mnsec ) mnsec = las;
        if ( las > 0 ) las = 0;
    }
}

inline LL solve( const Point& P ) {
    if ( !P.x ) return P.y > 0 ? -P.y * mnsec : -P.y * mxsec;
    if ( cvxA.size() == 1 ) return P ^ cvxA[0];
    Vector kp( P.x, P.y );
    if ( kp.x < 0 ) kp.x = -kp.x, kp.y = -kp.y;

    int s = int( cvxA.size() );
    if ( P.x > 0 ) { // up convex.
        int l = rpos, r = s;
        while ( l < r ) {
            int mid = l + r >> 1;
            if ( ( ( cvxA[mid] - cvxA[( mid + 1 ) % s] ) ^ kp ) > 0 )
                l = mid + 1;
            else r = mid;
        }
        l %= s;
        return P.x * cvxA[l].y - P.y * cvxA[l].x;
    } else { // down convex.
        int l = 0, r = rpos;
        while ( l < r ) {
            int mid = l + r >> 1;
            if ( ( ( cvxA[( mid + 1 ) % s] - cvxA[mid] ) ^ kp ) > 0 ) 
                l = mid + 1;
            else r = mid;
        }
        l %= s;
        return P.x * cvxA[l].y - P.y * cvxA[l].x;
    }
}

int main() {
    freopen( "area.in", "r", stdin );
    freopen( "area.out", "w", stdout );

    n = rint(), q = rint();
    rep ( i, 1, n ) A[i].x = rint(), A[i].y = rint();

    init();

    for ( Point P; q--; ) {
        P.x = rint(), P.y = rint();
        wint( lmax( solve( P ), 0 ) ), putchar( '\n' );
    }
    return 0;
}