zoukankan      html  css  js  c++  java
  • 【ZJOI2014】力

    Link:
    Luogu https://www.luogu.com.cn/problem/P3338


    Description

    [forall i=1,2,cdots,n ]

    ,求:

    [E_i~=~sumlimits_{j = 1}^{i - 1} frac{q_j}{(i - j)^2}~-~sumlimits_{j = i + 1}^{n} frac{q_j}{(i - j)^2} ]

    其中

    [n~in~mathbb{N}~cap~[1,10^5], ]

    [forall j~in~mathbb{N}~cap~[1,n]~,~quad q_i~in~mathbb{R}~cap~(0,10^9) ]



    Solution

    怎么说呢,

    [sum_{j = 1}^{i - 1} frac{q_j}{(i - j)^2} ]

    算是个卷积吧
    那右边那个反正直觉就是可以变成卷积
    虽然但是,不难发现把 (i) 改成倒过来的马上就得到一个卷积了
    (比如说,原来的 (i=n) 就变成 (i=1)

    更加具体地,左边是
    首先设 (f(j)=q_j) 以及 (g(x)=x^{-2})

    [sumlimits_{j=1}^{i-1}f(j)g(i-j) ]

    右边是,

    [sumlimits_{t=1}^{n-i}frac{q_{i+t}}{t^2} ]

    不妨设

    [h(x)quad s.t.quad h(n-i-t)=q_{i+t} ]

    ……至于下标不够怎么处理?
    补0就可以啦()

    一度有过这样的神奇想法
    “为什么输出答案的时候,左边和右边两个卷积不用先求和呢?”
    ……没救了。卷积


    经典错误了
    就是那个
    lim
    。。因为卷积所以至少不小于 (2n)
    然后因为
    for(0; <lim; )
    所以实际上 lim 不能比 (2n+1) 小。



    Code

    #include<cstdio>
    #include<algorithm>
    #include<cstdlib>
    #include<cstring>
    #include<iostream>
    #include<cctype>
    #include<cmath>
    
    using namespace std;
    
    const int MAXN = 4e5 + 10;
    const double tpi = acos(-1.0);
    
    int n, lim = 1, loglim = 0;
    int rev[MAXN];
    
    struct complex
    {
        double real, imag;
    
        complex(const double &aug1 = 0, const double &aug2 = 0)
        {
            real = aug1;
            imag = aug2;
        }
        
        complex operator * (complex b)
        {
            return complex(real * b.real - imag * b.imag, real * b.imag + imag * b.real);
        }
    
        complex operator + (complex b)
        {
            return complex(real + b.real, imag + b.imag);
        }
        
        complex operator - (complex b)
        {
            return complex(real - b.real, imag - b.imag);
        }
    
        void operator *= (complex b)
        {
            *this = *this * b;
        }
    
        void operator += (complex b)
        {
            real += b.real;
            imag += b.imag;
        }
    
        void operator -= (complex b)
        {
            real -= b.real;
            imag -= b.imag;
        }
    
        void operator *= (const double &x)
        {
            real *= x;
            imag *= x;
        }
    
    } f[MAXN], g[MAXN], h[MAXN], Wn[MAXN][2];
    
    void fft(complex *f, int opt)
    {
        static complex t;
    
        for (int i = 0; i < lim; ++i)
        {
            if (rev[i] > i) swap(f[i], f[rev[i]]);
        }
    
        for (int n, stats, half_n = 1; half_n < lim; half_n = n)
        {
            n = half_n << 1;
            stats = lim / n;
            for (int pos = 0; pos < lim; pos += n)
            {
                for (int k = 0; k < half_n; ++k)
                {
                    t = Wn[stats * k][opt] * f[pos + half_n + k];
                    f[pos + half_n + k] = f[pos + k] - t;
                    f[pos + k] += t;
                }
            }
        }
    
        if (opt == 1)
        {
            for (int i = 0; i < lim; ++i)
                f[i].real /= lim;
        }
    }
    
    void init()
    {
        scanf("%d", &n);
    
        for (int i = 1; i <= n; ++i)
        {
            scanf("%lf", &f[i].real);
            h[n - i].real = f[i].real;
            g[i].real = 1.0 / i / i;
        }
        
        while (lim <= (n << 1))
        {
            lim <<= 1;
            ++loglim;
        }
    
        for (int i = 0; i < lim; ++i)
        {
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (loglim - 1));
        }
    
        const double modpi = 2.0 * tpi / lim;
    
        double theta;
        for (int i = 0; i < lim; ++i)
        {
            theta = modpi * i;
            Wn[i][0].real = cos(theta);
            Wn[i][1].real = +Wn[i][0].real;
            Wn[i][0].imag = sin(theta);
            Wn[i][1].imag = -Wn[i][0].imag;
        }
    }
    
    void work()
    {
        fft(f, 0);
        fft(g, 0);
        fft(h, 0);
    
        for (int i = 0; i < lim; ++i)
        {
            f[i] *= g[i];
            h[i] *= g[i];
        }
    
        fft(f, 1);
        fft(h, 1);
    }
    
    int main()
    {
        init();
    
        work();
    
        for (int i = 1; i <= n; ++i)
        {
            printf("%.3f
    ", f[i].real - h[n - i].real);
        }
    
        system("pause");
    
        return 0;
    }
    
  • 相关阅读:
    Ubuntu安装adobe的Source Code Pro
    Oracle实现主键自增的几种方式
    Oracle主键自增
    activity 根据流程实例ID删除流程实例、删除流程部署
    解决报错:错误1130- Host xxx is not allowed to connect to this MariaDb server
    Idea-每次修改JS文件都需要重启Idea才能生效解决方法
    html中<a href> </a>的用法
    Mysql8_数据库基础操作
    java反射及Method的Invoke方法(转载)
    oracle查询数据插入时间
  • 原文地址:https://www.cnblogs.com/ccryolitecc/p/14002245.html
Copyright © 2011-2022 走看看