Bzoj3527 [Zjoi2014]力

Time Limit: 30 Sec  Memory Limit: 256 MBSec  Special Judge
Submit: 1841  Solved: 1095

Description

给出n个数qi，给出Fj的定义如下：

令Ei=Fi/qi，求Ei.

 

Input

第一行一个整数n。

接下来n行每行输入一个数，第i行表示qi。

n≤100000，0<qi<1000000000

 

 

Output

 n行，第i行输出Ei。与标准答案误差不超过1e-2即可。

Sample Input

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

Sample Output

-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

HINT

 

Source

数学 FFT

 

式子里的q[j]可以直接约掉。

前后两半式子看上去很相似，所以拆开考虑。

对于前半边：

　　分母上的(j-i)是等差递减的，离j越近就越小，分子上的qi下标是等差递增的，离j越近就越大。

　　↑多么像一个优美的卷积

　　然后构造出一个新多项式B，第i项的系数为1/(i+1)^2，和q[]数组做卷积，结果的第i位就是E[i+1]的前一半 （i下标从0开始）；接着把q数组倒序，多项式B整体取负，做出来的结果E[i]就是E[n-(i+1)-1]的后一半。

　　↑数学多么优美

　　写完之后的半个小时里各种调试，就是过不了样例……然后想起来【第i位就是E[i+1]的前一半】，啊，需要把答案整体右移一位。最左边的E[1]当然没有前一半可以加，最右边的E[n]当前没有后一半可以减，所以这两位置0

　　然后就妥妥地过了

　　之后我在想，这个右移多麻烦啊。

　　确实呢，多项式B中，只要把第i项的系数设为1/i^2(下标从0开始，所以第0项为0)，结果的第i位就是E[i]的结果咯

/*by SilverN*/
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<vector>
using namespace std;
const double pi=acos(-1.0);
const int mxn=800010;
struct com{
    double x,y;
    com operator + (com b){return (com){x+b.x,y+b.y};}
    com operator - (com b){return (com){x-b.x,y-b.y};}
    com operator * (com b){return (com){x*b.x-y*b.y,x*b.y+y*b.x};}
}a[mxn],b[mxn],c[mxn];
int n,l;
int rev[mxn];
void FFT(com *a,int flag){
    int i,j,k;
    for(i=0;i<n;i++)
        if(rev[i]>i)swap(a[rev[i]],a[i]);
    for(i=1;i<n;i<<=1){
        com wn=(com){cos(pi/i),flag*sin(pi/i)};
        for(j=0;j<n;j+=(i<<1)){
            com w=(com){1,0};
            for(k=0;k<i;k++,w=w*wn){
                com x=a[j+k],y=w*a[i+j+k];
                a[j+k]=x+y;
                a[i+j+k]=x-y;
            }
        }
    }
    if(flag==-1)for(i=0;i<n;i++) a[i].x/=n;
    return;
}
int len;
double q[mxn];
int main(){
    int i,j;
    scanf("%d",&len);
    for(i=0;i<len;i++)scanf("%lf",&q[i]);
    //
    int m=len<<1;
    for(n=1;n<m;n<<=1)l++;
    for(i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    //
    for(i=0;i<len;i++)a[i].x=q[i];
    for(i=1;i<len;i++)b[i].x=1.0/(double)(i)/(double)(i);
    //
    FFT(a,1);FFT(b,1);
    for(i=0;i<n;i++)c[i]=a[i]*b[i];
    //
    FFT(c,-1);
    
    
    memset(a,0,sizeof a);
    memset(b,0,sizeof b);
    for(i=0;i<len;i++)a[i].x=q[len-i-1];
    for(i=1;i<len;i++)b[i].x=-1.0/(double)(i)/(double)(i);
    
    FFT(a,1);FFT(b,1);
    for(i=0;i<n;i++)a[i]=a[i]*b[i];
//    FFT(c,-1);
    FFT(a,-1);
/*    
    for(i=n-1;i;i--){
        c[i]=c[i-1];
        a[i]=a[i-1];
    }
    c[0]=a[0]=(com){0,0};
*/
    for(i=0;i<len;i++)c[i]=c[i]+a[len-i-1];
    for(i=0;i<len;i++){
        printf("%.3f\n",c[i].x);
    }
    return 0;
}