bzoj3527: [Zjoi2014]力

我有优越感，我只跑了一次卷积~

推了一中午柿子：

其实就是求Ej=sigema(i<j)i qi/(i-j)^2 - sigema(j<i)i qi/(i-j)^2;

设n=5
E1= q1/(0^2) -q2/(1^2) -q3/(2^2) -q4/(3^2) -q5/(4^2)
E2= q1/(1^2) +q2/(0^2) -q3/(1^2) -q4/(2^2) -q5/(3^2)
E3= q1/(2^2) +q2/(1^2) +q3/(0^2) -q4/(1^2) -q5/(2^2)
E4= q1/(3^2) +q2/(2^2) +q3/(1^2) +q4/(0^2) -q5/(1^2)
E5= q1/(4^2) +q2/(3^2) +q3/(2^2) +q4/(1^2) +q5/(0^2)

A1~n=q1~qn B1~n=-1/((n-1)^2)~1/(0^2) B5视为0
实际上式子变成了这样
Ej=sigema(i<j)i -Ai*B(n-(j-i+1)+1) + sigema(j<i)i Ai*B(n-(i-j+1)+1);
Ej=sigema(i<j)i -Ai*B(n-j+i) + sigema(j<i)i Ai*B(n-i+j);

栗子
E1= +A1*B5 +A2*B4 +A3*B3 +A4*B2 +A5*B1
E2= -A1*B4 +A2*B5 +A3*B4 +A4*B3 +A5*B2
E3= -A1*B3 -A2*B4 +A3*B5 +A4*B4 +A5*B3
E4= -A1*B2 -A2*B3 -A3*B4 +A4*B5 +A5*B4
E5= -A1*B1 -A2*B2 -A3*B3 -A4*B4 +A5*B5

视为C6~C10 目前式子变成了这样
Cj=sigema(i<j-n)i -Ai*B(n-(j-n)+i) + sigema(j-n<i)i Ai*B(n-i+(j-n));
Cj=sigema(i<j-n)i -Ai*B(2*n-j+i) + sigema(j-n<i)i Ai*B(j-i);

令Bn~n*2=1/(0^2)~1/((n-1)^2)
那么-B(i)=B((2*n-1)-i+1)=B(2*n-i)
那 -B(2*n-j+i)=B(2*n-(2*n-j+i))=B(j-i)

Cj=sigema(i<j-n)i Ai*B(j-i) + sigema(j-n<i)i Ai*B(j-i);
C(j)=sigema(1~j-n)i A(i)*B(j-i)

C6= A1*B5 +A2*B4 +A3*B3 +A4*B2 +A5*B1
C7= A1*B6 +A2*B5 +A3*B4 +A4*B3 +A5*B2
C8= A1*B7 +A2*B6 +A3*B5 +A4*B4 +A5*B3
C9= A1*B8 +A2*B7 +A3*B6 +A4*B5 +A5*B4
C10= A1*B9 +A2*B8 +A3*B7 +A4*B6 +A5*B5

然后就FFT咯，结果模版打错了一次囧

update:我是神吧我写的是啥啊。。。我自己都不懂只会两次卷积做啊。。。。

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
using namespace std;
const double pi=acos(-1.0);

struct Complex
{
    double r,i;
    Complex(){}
    Complex(double _r,double _i){r=_r, i=_i;}
    friend Complex operator +(Complex x,Complex y){return Complex(x.r+y.r,x.i+y.i);}
    friend Complex operator -(Complex x,Complex y){return Complex(x.r-y.r,x.i-y.i);}
    friend Complex operator *(Complex x,Complex y){return Complex(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);}    
}A[810000],B[810000],C[810000];

int R[810000];
void fft(Complex *a,int n,int op)
{
    for(int i=0;i<n;i++)
        if(i<R[i])swap(a[i],a[R[i]]);
        
    for(int i=1;i<n;i*=2)
    {
        Complex wn(cos(pi/i),sin(pi*op/i));
        for(int j=0;j<n;j+=(i<<1))
        {
            Complex w(1,0);
            for(int k=0;k<i;k++,w=w*wn)
            {
                Complex a1=a[j+k],a2=a[j+k+i];
                a[j+k]  =a1+w*a2;
                a[j+k+i]=a1-w*a2;
            }
        }
    }
}

int n,m;
void yu()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%lf",&A[i].r);
    m=n*2-1;
    for(int i=1;i<=n;i++)
    {
        if(i==n)B[i].r=0;
        else    B[i].r=-1.0/(double(n-i)*double(n-i));
        if(i==1)B[n+i-1].r=0;
        else     B[n+i-1].r=1.0/(double(i-1)*double(i-1));
    }
    
    for(int i=1;i<=n;i++)A[i-1].r=A[i].r;
    for(int i=1;i<=m;i++)B[i-1].r=B[i].r;
    A[n].r=0;n--;
    B[m].r=0;m--;
}
int main()
{
    yu();
    int L=0,tt=n;
    m+=n;for(n=1;n<=m;n*=2)L++;
    for(int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    

    fft(A,n,1);fft(B,n,1);
    for(int i=0;i<n;i++)C[i]=A[i]*B[i];
    
    fft(C,n,-1);
    for(int i=tt;i<=tt*2;i++)printf("%lf
",C[i].r/double(n));
    return 0;
}

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
using namespace std;
const double pi=acos(-1.0);

struct complex
{
    double r,i;
    complex(){}
    complex(double R,double I){r=R,i=I;}
    friend complex operator +(complex x,complex y){return complex(x.r+y.r,x.i+y.i);}
    friend complex operator -(complex x,complex y){return complex(x.r-y.r,x.i-y.i);}
    friend complex operator *(complex x,complex y){return complex(x.r*y.r-x.i*y.i,x.r*y.i+y.r*x.i);}
}A[410000],B[410000],C[410000];
int Re[410000];
void fft(complex *a,int n,int op)
{
    for(int i=0;i<n;i++)
        if(i<Re[i])swap(a[i],a[Re[i]]);
    
    for(int i=1;i<n;i*=2)
    {
        complex wn(cos(pi/i),sin(op*pi/i));
        for(int j=0;j<n;j+=(i<<1))
        {
            complex w(1,0);
            for(int k=0;k<i;k++,w=w*wn)
            {
                complex t1=a[j+k],t2=a[j+k+i]*w;
                a[j+k]=t1+t2;
                a[j+k+i]=t1-t2;
            }
        }
    }
}

int n,pn; double q[110000];
void solve()
{
    memset(A,0,sizeof(A));
    memset(B,0,sizeof(B));
    memset(C,0,sizeof(C));
    A[0].r=q[0];
    for(int i=1;i<=pn;i++)
        A[i].r=q[i], B[i].r=1.0/(double(i)*double(i));
        
    fft(A,n,1),fft(B,n,1);
    for(int i=0;i<=n;i++)C[i]=A[i]*B[i];
    fft(C,n,-1);
    for(int i=0;i<=n;i++)C[i].r=C[i].r/double(n);
}
double as[110000];
int main()
{
    freopen("a.in","r",stdin);
    freopen("a.out","w",stdout);
    int m,L=0;
    scanf("%d",&n);n--,pn=n;
    for(int i=0;i<=n;i++)scanf("%lf",&q[i]);
    m=n*2;for(n=1;n<=m;n*=2)L++;
    for(int i=0;i<n;i++)Re[i]=(Re[i>>1]>>1)|((i&1)<<(L-1));
    
    memset(as,0,sizeof(as));
    solve();
    for(int i=0;i<=pn;i++)as[i]+=C[i].r;
    for(int i=0;i<=pn/2;i++)swap(q[i],q[pn-i]);
    solve();
    for(int i=0;i<=pn;i++)as[i]-=C[pn-i].r;
    
    for(int i=0;i<=pn;i++)printf("%.5lf
",as[i]);
    return 0;
}