HDU
来源:互联网 发布:javascript长度函数 编辑:程序博客网 时间:2024/06/03 20:27
题目链接点这里
过程看这里http://blog.csdn.net/u013368721/article/details/53001532
,,然而,,这位大佬的程序也gg了,,可能改过数据了,,不过思路是这样的,,
然而我与标程对拍,,答案总是差0.00几。。反正,,我写的fft自带低精度,大内存,高常数。。,,以后再看吧
#include<iostream>#include<cstdio>#include<math.h>#include<algorithm>#include<map>#include<set>#include<bitset>#include<stack>#include<queue>#include<string.h>#include<cstring>#include<vector>#include<time.h>#include<stdlib.h>using namespace std;#define INF 0x3f3f3f3f#define INFLL 0x3f3f3f3f3f3f3f3f#define FIN freopen("input.txt","r",stdin)#define mem(x,y) memset(x,y,sizeof(x))typedef unsigned long long ULL;typedef long long LL;#define fuck(x) cout<<"x"<<endl;#define lson l,m,rt<<1#define rson m+1,r,rt<<1|1typedef pair<pair<int,int>,int> PIII;typedef pair<int,int> PII;const double eps=1e-100;const int MX=1111111;int P;int G;double a[MX];//fftconst double pi = acos(-1.0);int len,res[MX],mx;//开大4倍struct Complex{ double r,i; Complex(double r=0,double i=0):r(r),i(i) {}; Complex operator+(const Complex &rhs) { return Complex(r + rhs.r,i + rhs.i); } Complex operator-(const Complex &rhs) { return Complex(r - rhs.r,i - rhs.i); } Complex operator*(const Complex &rhs) { return Complex(r*rhs.r - i*rhs.i,i*rhs.r + r*rhs.i); }} va[MX],vb[MX];void rader(Complex F[],int len) //len = 2^M,reverse F[i] with F[j] j为i二进制反转{ int j = len >> 1; for(int i = 1; i < len - 1; ++i) { if(i < j) swap(F[i],F[j]); // reverse int k = len>>1; while(j>=k) { j -= k; k >>= 1; } if(j < k) j += k; }}void FFT(Complex F[],int len,int t){ rader(F,len); for(int h=2; h<=len; h<<=1) { Complex wn(cos(-t*2*pi/h),sin(-t*2*pi/h)); for(int j=0; j<len; j+=h) { Complex E(1,0); //旋转因子 for(int k=j; k<j+h/2; ++k) { Complex u = F[k]; Complex v = E*F[k+h/2]; F[k] = u+v; F[k+h/2] = u-v; E=E*wn; } } } if(t==-1) //IDFT for(int i=0; i<len; ++i) F[i].r/=len;}/*inline void FFT(Complex *a,int n,int r){ rader(a,len); for (int i=1; i<n; i<<=1) { Complex wn(cos(pi/i),r*sin(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 x=a[j+k],y=w*a[j+k+i]; a[j+k]=x+y; a[j+k+i]=x-y; } } } if (r==-1) for (int i=0; i<n; i++) a[i].r/=n;}*/void Conv(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);}int change[MX];double ans[MX];double ss(int ret){ return pow(2,sin(2*pi*ret/P)*sin(2*pi*ret/P)*sin(2*pi*ret/P));}void gao(){ mx=P-2+P-2; len=1; while(len<mx)len<<=1;//mx为卷积后最大下标 int ret=1; va[P-2-0].r=a[1]; vb[0].r=ss(ret); change[0]=1; for(int i=1; i<P-1; i++) { ret=(ret*G)%P; change[i]=ret; va[P-2-i].r=a[ret]; va[P-2-i].i=0; vb[i].r=ss(ret); vb[i].i=0; } for(int i=P-1; i<len; i++)va[i].i=va[i].r=vb[i].i=vb[i].r=0; Conv(va,vb,len); ans[1]=va[P-2].r; ans[0]=a[P-1]; for(int i=0; i<P-2; i++) { ans[change[i+1]]=va[i].r+va[i+P-1].r; ans[0]+=a[i+1]; } for(int i=0; i<=P-1; i++) { printf("%.3f ",ans[i]+a[0]); } puts("");}int main(){ FIN; freopen("output1.txt","w",stdout); while(~scanf("%d",&P)) { if(P==13)G=2; else if(P==103)G=5; else G=2; for(int i=0; i<P; i++)scanf("%lf",&a[i]); gao(); } return 0;}
阅读全文
0 0
- hdu
- hdu
- HDU
- hdu ()
- hdu
- hdu
- HDU
- HDU
- hdu
- hdu
- HDU
- Hdu
- hdu
- hdu-
- hdu
- hdu
- hdu
- HDU
- centos7 编译安装tensorflow1.2
- jpa nulls last 不起作用
- Java注解(Annotation)详解(四)——注解反射生成SQL语句
- 564. Find the Closest Palindrome
- 省市县三级自身一对一关联的表,联动查询SQL
- HDU
- 移动开发----biu,biu,一个有趣的EditText
- Android源码-高质量开发库
- 【c语言】打开一个客户端socket描述符
- 包装类和正则表达式
- 图像增强(直方图均衡化、拉普拉斯、Log、Gamma)
- 控制开关的printf
- zoj1149 Dividing 动态规划
- RB_tree的插入