[Codeforces g100016B]Command Post (中文：圈地問題)

題目敘述
傳code

建議先將所有觀察站位置複製，變成$2N$個觀察站
處理環狀模型會方便許多

由於最後圍籬會形成一個多邊形，我們先枚舉觀察站位置，並強迫在這個地方蓋觀察站
因此這個位置就一定是多邊形的一個頂點，這樣剩下來的觀察站就有特定順序了 (問題變成非環形)，依序給它們編號$1\sim N-1$ (也就是，假設枚舉到第$i$個觀察站，那麼第$j$個觀察站的編號就是$j-i$)

說明：
$G[i]$為編號$i$的觀察站位置所在的角度 (參見原題定義)

先來考慮一個很樸素的$DP$
$DP[n][i]$代表在觀察站位置$1\sim i$選$n$個蓋觀察站 (有個限制是觀察站$i$一定要蓋)，對圓心而言所能覆蓋的最大面積 (就是兩兩點間的弦+它們的兩條半徑形成的很多三角形的面積的總和)
因此我們就有以下遞推式：
$DP[n][i]=\max_{j=0}^{i-1}(DP[n-1][j]+\frac{\sin(G[i]-G[j])}{2},G[i]-G[j]\leq\pi)$
邊界條件：
$DP[0][0]=0$

乍看之下很像斜率優化$DP$？
其實不然
因為它無法表示成一條條的直線QAQ

不過這個$DP$有一個性質
令$j_{i}$為第一個$j$使得$DP[n-1][j]+\frac{\sin(G[i]-G[j])}{2}$最大，$G[i]-G[j]\leq\pi$
則我們有$j_{i}\leq j_{i+1}$
限於篇幅，性質證明請看這裡

因此當我們算出$j_{k}$時，同時間也知道了$j_{i}\leq j_{k},\forall i\leq k$和$j_{i}\geq j_{k},\forall i\geq k$

利用這個性質，我們可以按照$n$慢慢遞增的順序計算$DP$
然後對於$DP[l\sim r]$，假設目前知道了$l_{b}\leq j_{l\sim r}\leq r_{b}$
那麼我們可以取$mid=\left\lfloor\frac{l+r}{2}\right\rfloor$，暴力枚舉$l_{b}\sim r_{b}$，算出$j_{mid}$ (因為我們已經知道$l_{b}\leq j_{mid}\leq r_{b}$)
於是
對於$DP[l\sim mid-1]$，我們知道了$l_{b}\leq j_{l\sim mid-1}\leq j_{mid}$
對於$DP[mid+1\sim r]$，我們知道了$j_{mid}\leq j_{mid+1\sim r}\leq r_{b}$

將這兩個子問題遞迴計算即可

可以證明，整體時間複雜度為$O(N\log N)$ (證明略去，但還是建議讀者自行把code打出來並證明一下，不會很難的)
因此算到$n=K$需要$O(KN\log N)$的時間

所以我們得到了$O(KN^{2}\log N)$的算法

當然，這樣還不夠快

可以發現，我們的計算方向是$n$遞增
因此與其每枚舉一個觀察站位置就花$O(KN\log N)$的時間去計算
不如全部的枚舉情況一起算(?)

令$A[n][i][j]$為觀察站$i+1\sim j$中選$n$個時，最大面積是多少 (觀察站$i$和$j$是強迫要選的)
因此，我們有了$A[a+b][i][j]=\max_{k=i+1}^{j-1}(A[a][i][k]+A[b][k][j])$
注意到，這很像矩陣乘法
只是這裡將元素的相乘變成相加，相加變成取最大值
而且利用前述性質的拓展 (請讀者自行思考怎麼拓展性質，不會很困難的)，這個特殊的矩陣乘法可以做到$O(N^{2}\log N)$！
或許我們可以利用快速冪(?)

我們可以先算出$A[1]$
利用$A[1]*A[1]$我們可以算出$A[2]$
利用$A[2]*A[2]$我們可以算出$A[4]$
利用$A[4]*A[4]$我們可以算出$A[8]$
......
任何$2$的冪次方都可以藉由$O(\log X)$次矩陣乘法算出來！

非$2$的冪次方怎麼辦呢？
以$19$次方為例
我們可以先用$A[1]$和$A[2]$相乘算出$A[3]$
然後再用$A[3]$和$A[16]$算出$A[19]$！
一般的
把數字轉換成二進位 (例如$19$就是$10011$)
右邊數來第$x$位是$1$就把矩陣乘上$A[2^{x-1}]$
任何冪次方都可以藉由$O(\log X)$次矩陣乘法算出來！

呵呵，我們講到哪裡了？
對啊，這樣把$A[K]$算出來
然後枚舉$i$找出最大的$A[K][i][i+N]$ (因為要完整環繞一圈變成多邊形)
就知道最大面積多少了！

可是題目要我們輸出方案
簡單！(拍桌)
記錄下中間結果再逆推回去就好啦！

時間複雜度：$O(N^{2}\log N\log K)$

p.s.你以為這麼容易嗎？XD
p.s.這份code會MLE (還沒結束?!)，我們必須有效利用未使用的儲存空間，將記憶體用量減為$\frac{1}{4}$

code:

#include<cstdio>
#include<cassert>
#include<vector>
#include<algorithm>
#include<cmath>
using namespace std;
const double INF=1e100,EPS=1e-7,PI=acos(-1.0);
bool getmax(double &a,const double &b){if(b>a){a=b;return true;}return false;}
bool getmax(double &a,const double &b1,const double &b2)
{
    if(b1==-INF||b2==-INF)return false;
    if(b1+b2>a){a=b1+b2;return true;}
    return false;
}
int N,K;
struct Matrix
{
    double S[2000][2000];
    int MID[2000][2000];
    Matrix *ch[2];
    Matrix(Matrix *a,Matrix *b):ch{a,b}{}
};
//Matrix *Mul(Matrix *a,Matrix *b)
//{
//    Matrix *ans=new Matrix(a,b);
//    for(int i=0;i<N*2;i++)for(int j=0;j<N*2;j++)
//    {
//        double &v=ans->S[i][j]=-INF;
//        int &mid=ans->MID[i][j]=-1;
//        for(int k=0;k<N*2;k++)if(getmax(v,a->S[i][k],b->S[k][j]))mid=k;
//    }
//    return ans;
//}
void Solve(const int from,const int l,const int r,const int bl,const int br,Matrix *ans,Matrix *a,Matrix *b)
{
    if(l>r)return;
    const int mid=(l+r)/2;
    int &bm=ans->MID[from][mid]=bl;
    double &v=ans->S[from][mid]=-INF;
    for(int i=bl;i<=br&&i<=mid;i++)
    {
        if(getmax(v,a->S[from][i],b->S[i][mid]))bm=i;
    }
    Solve(from,l,mid-1,bl,bm,ans,a,b),Solve(from,mid+1,r,bm,br,ans,a,b);
}
Matrix *Mul(Matrix *a,Matrix *b)
{
    Matrix *ans=new Matrix(a,b);
    for(int i=0;i<N*2;i++)
    {
        for(int j=0;j<i;j++)ans->S[i][j]=-INF;
        Solve(i,i,N*2-1,i,N*2-1,ans,a,b);
    }
    return ans;
}
void ExtractPath(Matrix *o,const int l,const int r,vector<int>&output)
{
    if(!o->ch[0])return;
    const int mid=o->MID[l][r];
    assert(mid!=-1);
    ExtractPath(o->ch[0],l,mid,output);
    output.push_back(mid);
    ExtractPath(o->ch[1],mid,r,output);
}
struct Opl
{
    double ANGLE;
    int ID;
    bool operator<(const Opl &o)const{return ANGLE<o.ANGLE;}
};
Opl OPL[2000];
double DP[2000][2000];
int main()
{
    freopen("post.in","r",stdin);
//    freopen("in.txt","r",stdin);
    while(scanf("%d%d",&N,&K)==2)
    {
        for(int i=0;i<N;i++)scanf("%lf",&OPL[i].ANGLE),OPL[i].ID=i+1;
        sort(OPL,OPL+N);
        for(int i=0;i<N;i++)OPL[N+i].ID=OPL[i].ID,OPL[N+i].ANGLE=OPL[i].ANGLE+2.0*PI;
        static Matrix *m[12];
        vector<Matrix*>trash;
        trash.push_back(m[0]=new Matrix(NULL,NULL));
        for(int i=0;i<N*2;i++)for(int j=0;j<N*2;j++)
        {
            double &v=m[0]->S[i][j]=-INF;
            if(i<j&&OPL[j].ANGLE-OPL[i].ANGLE-EPS<=PI)v=0.5*sin(OPL[j].ANGLE-OPL[i].ANGLE);
        }
        int bound=0;
        for(int i=1;(1<<i)<=K;i++)trash.push_back(m[i]=Mul(m[i-1],m[i-1])),bound=i;
        Matrix *mans=NULL;
        for(int i=bound;i>=0;i--)if(K&(1<<i))
        {
            if(!mans)mans=m[i];
            else trash.push_back(mans=Mul(mans,m[i]));
        }
        double ans=0.0;
        int best_loc=-1;
        for(int i=0;i<N;i++)if(getmax(ans,mans->S[i][i+N]))best_loc=i;
//        printf("%.3lf\n",ans);
        if(best_loc==-1)puts("0");
        else
        {
            vector<int>path;
            path.push_back(best_loc);
            ExtractPath(mans,best_loc,best_loc+N,path);
            for(int &v:path)v%=N;
            sort(path.begin(),path.end());
            for(int i=0;i<(int)path.size();i++)
            {
                if(i)putchar(' ');
                printf("%d",OPL[path[i]].ID);
            }
            puts("");
            assert((int)path.size()==K);
        }
        for(Matrix *v:trash)delete v;
    }
    return 0;
}