RMQ的ST演算法

RMQ 即 range maximum/minimum query

　　1、樸素（即搜尋），O(n)-O(qn) online。

　　2、線段樹，O(n)-O(qlogn) online。

　　3、ST（實質是動態規劃），O(nlogn)-O(1) online。

　　ST演算法（Sparse Table），以求最大值為例，設d[i,j]表示[i,i+2^j-1]這個區間內的最大值，那麼在詢問到[a,b]區間的最大值時答案就是max(d[a,k], d[b-2^k+1,k])，其中k是滿足2^k<=b-a+1(即長度)的最大的k,即k=[ln(b-a+1)/ln(2)]。

　　d的求法可以用動態規劃，d[i,j]=max(d[i,j-1],d[i+2^(j-1),j-1])。

　　4、RMQ標準演算法：先規約成LCA（Lowest Common Ancestor），再規約成約束RMQ，O(n)-O(1) online。

 

　　首先根據原數列，建立笛卡爾樹，從而將問題線上性時間內規約為LCA問題。LCA問題可以線上性時間內規約為約束RMQ，也就是數列中任意兩個相鄰的數的差都是+1或-1的RMQ問題。約束RMQ有O(n)-O(1)的線上解法，故整個演算法的時間複雜度為O(n)-O(1)。

 

 ST演算法　　來看一下ST演算法是怎麼實現的（以最大值為例）：

　　首先是預先處理，用一個DP解決。設a是要求區間最值的數列，f[i,j]表示從第i個數起連續2^j個數中的最大值。例如數列3 2 4 5 6 8 1 2 9 7 ,f[1,0]表示第1個數起，長度為2^0=1的最大值，其實就是3這個數。f[1,2]=5，f[1,3]=8，f[2,0]=2，f[2,1]=4……從這裡可以看出f[i,0]其實就等於a。這樣，DP的狀態、初值都已經有了，剩下的就是狀態轉移方程。我們把f[i,j]（j≥1）平均分成兩段（因為j≥1時，f[i,j]一定是偶數個數字），從i到i+2^(j-1)-1為一段，i+2^(j-1)到i+2^j-1為一段（長度都為2^（j-1））。用上例說明，當i=1，j=3時就是3,2,4,5
和6,8,1,2這兩段。f就是這兩段的最大值中的最大值。於是我們得到了動規方程F[i,j]=max(F[i,j-1],F[i+2^(j-1),j-1])。

　　接下來是得出最值，也許你想不到計算出f有什麼用處，一般毛想想計算max還是要O(logn)，甚至O(n)。但有一個很好的辦法，做到了O(1)。還是分開來。如在上例中我們要求區間[2,8]的最大值，就要把它分成[2,5]和[5,8]兩個區間，因為這兩個區間的最大值我們可以直接由f[2,2]和f[5,2]得到。擴充到一般情況，就是把區間[l,r]分成兩個長度為2^n的區間（保證有f對應）。直接給出運算式：

　　k:=trunc(ln(r-l+1)/ln(2));

　　ans:=max(F[l,k],F[r-2^k+1,k]);

　　這樣就計算了從l開始，長度為2^k的區間和從r-2^k+1開始長度為2^k的區間的最大值（運算式比較煩瑣，細節問題如加1減1需要仔細考慮），二者中的較大者就是整個區間[l,r]上的最大值。

應用：

JOJ 2660 河床

時間上明顯沒有直接單調隊列維護最大最小快，可能是運算log上浪費了時間。

對log最佳化了下，直接位元運算取ln2稍微快了些

#include <cstdio>#include <cmath>#define max(a,b) (a>b?a:b)#define min(a,b) (a<b?a:b)const int maxn=30005;const double ln2=log(2.0);int dpmax[maxn][16],dpmin[maxn][16],w[maxn];int n;int getln2 (int x){    int i;    for(i=0 ; x ; ++i,x>>=1);    return i-1;}void rmqinit (){    for (int i=0 ; i<n ; ++i)    dpmax[i][0]=dpmin[i][0]=w[i];    int m=getln2(n);    for (int i=1 ; i<=m ; ++i)    {        for (int j=n-1 ; j>=0 ; --j)        {            dpmax[j][i]=dpmax[j][i-1];//F[i,j]=max(F[i,j-1],F[i+2^(j-1),j-1])。            if(j+(1<<(i-1))<=n)dpmax[j][i]=max(dpmax[j][i],dpmax[j+(1<<(i-1))][i-1]);            dpmin[j][i]=dpmin[j][i-1];            if(j+(1<<(i-1))<=n)dpmin[j][i]=min(dpmin[j][i],dpmin[j+(1<<(i-1))][i-1]);        }    }}int rmqmin(int l,int r){    int m=getln2(r-l+1);//k:=trunc(ln(r-l+1)/ln(2));    return min(dpmin[l][m],dpmin[r-(1<<m)+1][m]);//ans:=min(F[l,k],F[r-2^k+1,k]);}int rmqmax(int l,int r)//O(1){    int m=getln2(r-l+1);    return max(dpmax[l][m],dpmax[r-(1<<m)+1][m]);}int k;int main (){    while (~scanf("%d%d",&n,&k),(n||k))    {                //printf("%d\n",getln2(6));        for (int i=0 ; i<n ; ++i)        {            scanf("%d",w+i);        }        rmqinit();        int start=0,end=0;        int ans=0;        while(true)        {            int tmpmax=rmqmax(start,end);            int tmpmin=rmqmin(start,end);            //printf("%d  %d  start=%d e=%d\n",tmpmin,tmpmax,start,end);            while(tmpmax-tmpmin<=k && end<n)            {                end++;                tmpmax=rmqmax(start,end);                tmpmin=rmqmin(start,end);            }            //printf("%d  %d  start=%d e=%d\n",tmpmin,tmpmax,start,end);            ans=max(ans,end-start);            tmpmax=rmqmax(start,end);            tmpmin=rmqmin(start,end);            //printf("%d  %d  start=%d e=%d\n",tmpmin,tmpmax,start,end);            while(tmpmax-tmpmin>k && start <=end)            {                start++;                tmpmax=rmqmax(start,end);                tmpmin=rmqmin(start,end);            }            //printf("%d\n",ans);            if(start+ans>n || end>=n )break;        }        printf("%d\n",ans);    }    return 0;}