[jzoj 6087] [GDOI2019模拟2019.3.26] 获取名额 解题报告 (泰勒展开+RMQ+精度)

题目链接：

https://jzoj.net/senior/#main/show/6087

题目：

题解：

• 只需要统计$\prod_{i=l}^r (1-\frac{a_i}{x})$
• =$exp(\sum_{i=l}^r ln(1-\frac{a_i}{x}))(x>a_i)$
• 我们可以把$ln(1-x)|x<1|$泰勒展开，得到$-\sum_{i=1}^{∞}\frac{x^i}{i}=0-\frac{x}{1}-\frac{x^2}{2}-\frac{x^3}{3}-...$
• 那么里面化简得到$-\sum_{i=1}^{r} \sum_{j=1}^∞\frac{a_i^j}{x^j*j})=-\sum_{i=1}^∞\frac{1}{i*x^i}\sum_{j=l}^{r}a_j^i$
• 实测$∞$取$50$就可以，预处理前缀和来快速计算$\sum_{j=l}^{r}a_j^i$
• 注意把$x=x/max(a_i)$,$a_i=a_i/max(a_i)$，防止爆$double$的范围
• 但是发现精度会有很大问题
• 当$\frac{a_i}{x}$较大时，$ln$里面的东西会比较接近$0$，很难保证精度
• 为$\frac{a_i}{x}$设置一个$lim$，当$\frac{a_i}{x}>lim$时直接算，再递归左右区间
• $lim$取$0.5​$是可以过的
• 寻找区间内最大的$\frac{a_i}{x}$用$RMQ$
• 参考自大米饼大佬的博客

代码：

#include<bits/stdc++.h>
using namespace std;
typedef double db;

const int N=6e5+15;
const int M=50;
int n,q;
int bin[21],lg[N],f[N][21];
db res;
db a[N],s[M+1][N];
inline char gc()
{
  static char*p1,*p2,s[1000000];
  if(p1==p2)p2=(p1=s)+fread(s,1,1000000,stdin);
  return(p1==p2)?EOF:*p1++;
}
inline int read()
{
    char ch=gc();int s=0;
    while (ch<'0'||ch>'9') ch=gc();
    while (ch>='0'&&ch<='9') {s=(s<<3)+(s<<1)+ch-'0';ch=gc();}
    return s;
}
int Max(int x,int y) {return a[x]>a[y]?x:y;}
void pre()
{
    for (int i=bin[0]=1;i<=20;bin[i]=bin[i-1]<<1,i++);
    lg[0]=-1;
    for (int i=1;i<=n;i++) f[i][0]=i,lg[i]=lg[i>>1]+1;
    for (int i=1;i<=20;i++)
        for (int j=1;j+bin[i]-1<=n;j++)
            f[j][i]=Max(f[j][i-1],f[j+bin[i-1]][i-1]);
}
int ask(int l,int r)
{
    int t=lg[r-l+1];
    return Max(f[l][t],f[r-bin[t]+1][t]);
}
db x;
db cal(int l,int r)
{
    db re=0,t=1;
    for (int i=1;i<=M;i++)
    {
        t*=x;
        re-=(s[i][r]-s[i][l-1])/(i*t);
    }
    return exp(re);
}
void solve(int l,int r)
{
    int mid=ask(l,r);
    if (a[mid]/x<0.5) {res*=cal(l,r);return;}
    res*=1-a[mid]/x;
    if (l<mid) solve(l,mid-1);
    if (r>mid) solve(mid+1,r);
}
int main()
{
    freopen("orz.in","r",stdin);
    freopen("orz.out","w",stdout);
    n=read();q=read();
    db mx=0;
    for (int i=1;i<=n;i++) a[i]=read(),mx=max(mx,a[i]);
    for (int i=1;i<=n;i++) a[i]/=mx;
    pre();
    for (int i=1;i<=n;i++)
    {
        db t=1;
        for (int j=1;j<=M;j++)
        {
            t*=a[i];
            s[j][i]=s[j][i-1]+t;
        }
    }
    while (q--)
    {
        int l=read(),r=read();
        x=read()/mx;
        res=1;solve(l,r);
        printf("%.10f\n",1-res);
    }
    return 0;
}