您的位置:首页 > 其它

hdu 1402 (FFT 模版)

2013-08-26 13:03 267 查看
/*
algorithm : High-Precision FFT

Memory 7092K	Time 359MS	Language G++

code by : zhyu
*/
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#define N 200005
#define pi acos(-1.0) // PI值
using namespace std;
struct complex
{
double r,i;
complex(double real=0.0,double image=0.0){
r=real;	i=image;
}
// 以下为三种虚数运算的定义
complex operator + (const complex o){
return complex(r+o.r,i+o.i);
}
complex operator - (const complex o){
return complex(r-o.r,i-o.i);
}
complex operator * (const complex o){
return complex(r*o.r-i*o.i,r*o.i+i*o.r);
}
}x1
,x2
;
char a[N/2],b[N/2];
int sum
; // 结果存在sum里
void brc(complex *y,int l) // 二进制平摊反转置换 O(logn)
{
register int i,j,k;
for(i=1,j=l/2;i<l-1;i++)
{
if(i<j)	swap(y[i],y[j]); // 交换互为下标反转的元素
// i<j保证只交换一次
k=l/2;
while(j>=k) // 由最高位检索,遇1变0,遇0变1,跳出
{
j-=k;
k/=2;
}
if(j<k)	j+=k;
}
}
void fft(complex *y,int l,double on) // FFT O(nlogn)
// 其中on==1时为DFT,on==-1为IDFT
{
register int h,i,j,k;
complex u,t;
brc(y,l); // 调用反转置换
for(h=2;h<=l;h<<=1) // 控制层数
{
// 初始化单位复根
complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
for(j=0;j<l;j+=h) // 控制起始下标
{
complex w(1,0); // 初始化螺旋因子
for(k=j;k<j+h/2;k++) // 配对
{
u=y[k];
t=w*y[k+h/2];
y[k]=u+t;
y[k+h/2]=u-t;
w=w*wn; // 更新螺旋因子
} // 据说上面的操作叫蝴蝶操作…
}
}
if(on==-1)	for(i=0;i<l;i++)	y[i].r/=l; // IDFT
}
int main(void)
{
int l1,l2,l;
register int i;
while(scanf("%s%s",a,b)!=EOF)
{
l1=strlen(a);
l2=strlen(b);
l=1;
while(l<l1*2 || l<l2*2)	l<<=1; // 将次数界变成2^n
// 配合二分与反转置换
for(i=0;i<l1;i++) // 倒置存入
{
x1[i].r=a[l1-i-1]-'0';
x1[i].i=0.0;
}
for(;i<l;i++)	x1[i].r=x1[i].i=0.0;
// 将多余次数界初始化为0
for(i=0;i<l2;i++)
{
x2[i].r=b[l2-i-1]-'0';
x2[i].i=0.0;
}
for(;i<l;i++)	x2[i].r=x2[i].i=0.0;
fft(x1,l,1); // DFT(a)
fft(x2,l,1); // DFT(b)
for(i=0;i<l;i++)	x1[i]=x1[i]*x2[i]; // 点乘结果存入a
fft(x1,l,-1); // IDFT(a*b)
for(i=0;i<l;i++)	sum[i]=x1[i].r+0.5; // 四舍五入
for(i=0;i<l;i++) // 进位
{
sum[i+1]+=sum[i]/10;
sum[i]%=10;
}
l=l1+l2-1;
while(sum[l]<=0 && l>0)	l--; // 检索最高位
for(i=l;i>=0;i--)	putchar(sum[i]+'0'); // 倒序输出
putchar('\n');
}
return 0;
}


转自:http://zhyu.me/acm/hdu-1402.html

蝶形算法,Orz。。。

http://blog.sina.com.cn/s/blog_592c83c40100odna.html

二分平摊反转置换,

/article/4948487.html

Mycode:

hdu 4609 3-idiots

#include <cstdio>
#include <deque>
#include <set>
#include <string>
#include <map>
#include <vector>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cmath>
#include <queue>
using namespace std;
typedef long long LL;
const int maxn = 200005;
const double Pi = acos(-1.0);

struct complex{
double r,i;
complex(double _r=0, double _i=0){
r = _r;
i = _i;
}

complex operator +(const complex &b){
return complex(r+b.r,i+b.i);
}

complex operator -(const complex &b){
return complex(r-b.r,i-b.i);
}

complex operator *(const complex &b){
return complex(r*b.r-i*b.i,r*b.i+i*b.r);
}
};

static complex x[maxn*2],y[maxn*2];

void brc(complex y[],int len){
int i,j,k;
for(i=1,j=len/2;i<len-1;i++){
if(i<j)swap(y[i],y[j]);
k = len / 2;
while(j >= k){
j -= k;
k /= 2;
}
if(j < k)j += k;
}
}

void FFT(complex y[],int len,int on){
brc(y,len);
for(int h=2;h<=len;h<<=1){
complex wn(cos(-on*2*Pi/h),sin(-on*2*Pi/h));
for(int j=0;j<len;j+=h){
complex w(1,0);
for(int k=j;k<j+h/2;k++){
complex u = y[k];
complex t = w * y[k+h/2];
y[k] = u + t;
y[k+h/2] = u - t;
w = w * wn;
}
}
}
if(on == -1){
for(int i=0;i<len;i++){
y[i].r /= len;
}
}
}

/* a (*) b = IDFT(DFT(a).DFT(b)) */
void mul(complex a[],complex b[],int len){
FFT(a,len,1);
FFT(b,len,1);
for(int i=0;i<len;i++){
a[i] = a[i] * b[i];
}
FFT(a,len,-1);
}

void mul(LL a[],LL b[],LL c[],int l){
int len = 1;
while(len < l)len *= 2;
len *= 2;
for(int i=0;i<l;i++){
x[i] = complex(a[i],0);
y[i] = complex(b[i],0);
}
for(int i=l;i<len;i++)x[i] = y[i] = complex();
mul(x,y,len);
for(int i=0;i<=2*l;i++){
c[i] = (LL)(x[i].r+0.5);
}
}

void printarray(LL a[],int n){
for(int i=0;i<n;i++){
printf("%lld ",a[i]);
}
printf("\n");
}

int a[maxn];
LL num[maxn*2];
LL sum[maxn*2];
int main(){
int T;
scanf("%d",&T);
while(T --){
int n;
scanf("%d",&n);
memset(num,0,sizeof(num));
for(int i=0;i<n;i++){
scanf("%d",&a[i]);
num[(a[i])] ++;
}
sort(a,a+n);
int _max = a[n-1] + 1;

//printarray(num,_max);

mul(num,num,num,_max);
_max *= 2;

//printarray(num,_max+1);

for(int i=0;i<n;i++){
num[a[i]+a[i]] --;
}
for(int i=1;i<=_max;i++){
num[i] /= 2;
}

sum[0] = 0;
for(int i=1;i<=_max;i++){
sum[i] = sum[i-1] + num[i];
}

LL ans = 0;
for(int i=0;i<n;i++){
ans += sum[_max] - sum[a[i]];
ans -= (LL)(n-1-i)*i;
ans -= (n-1);
ans -= (LL)(n-1-i)*(n-2-i)/2;
}

LL tot = (LL)n*(n-1)*(n-2)/6;
printf("%.7lf\n",1.0*ans/tot);
}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: