我居然不会写多项式卷积了,把$a[i]=a[i]*b[i]$写成了$a[i].x=a[i].x*b[i].x$还调了好长时间,有毒
把式子写成$C^2(x)-A^2(x)-B^2(x)=0$的形式,然后牛顿迭代即可
1 #include<cmath> 2 #include<cstdio> 3 #include<cctype> 4 #include<cstring> 5 #include<algorithm> 6 using namespace std; 7 const int N=800005; 8 const double pai=acos(-1),inf=1e9,eps=1e-6; 9 struct cpx 10 { 11 double x,y; 12 }a[N],b[N],c[N]; 13 cpx operator + (cpx a,cpx b) 14 { 15 return (cpx){a.x+b.x,a.y+b.y}; 16 } 17 cpx operator - (cpx a,cpx b) 18 { 19 return (cpx){a.x-b.x,a.y-b.y}; 20 } 21 cpx operator * (cpx a,cpx b) 22 { 23 double x1=a.x,y1=a.y,x2=b.x,y2=b.y; 24 return (cpx){x1*x2-y1*y2,x1*y2+y1*x2}; 25 } 26 int la,lb,lc,n,m,nm,rev[N]; 27 double ll,rr,Sin[32],Cos[32],f[N],d[N]; 28 void Read(double &x) 29 { 30 int ret=0; x=0; 31 char ch=getchar(); 32 while(!isdigit(ch)) 33 ch=getchar(); 34 while(isdigit(ch)) 35 ret=(ret<<1)+(ret<<3)+(ch^48),ch=getchar(); 36 x=ret; return ; 37 } 38 void Prework() 39 { 40 m=2*(max(la,max(lb,lc))+1),n=1; while(n<=m) n<<=1; 41 for(int i=0;i<=n;i++) 42 rev[i]=(rev[i>>1]>>1)+(i&1)*(n>>1); 43 for(int i=1;i<=24;i++) 44 Sin[i]=sin(2*pai/(1<<i)),Cos[i]=cos(2*pai/(1<<i)); 45 } 46 void Trans(cpx *cop,int len,int typ) 47 { 48 register int i,j,k; 49 for(i=0;i<=len;i++) 50 if(rev[i]>i) swap(cop[rev[i]],cop[i]); 51 for(i=2;i<=len;i<<=1) 52 { 53 int lth=i>>1,lgg=log2(i); 54 cpx omg=(cpx){Cos[lgg],typ*Sin[lgg]}; 55 for(j=0;j<len;j+=i) 56 { 57 cpx ori={1,0}; 58 for(k=j;k<j+lth;k++,ori=ori*omg) 59 { 60 cpx tmp=ori*cop[k+lth]; 61 cop[k+lth]=cop[k]-tmp,cop[k]=cop[k]+tmp; 62 } 63 } 64 } 65 if(typ==-1) 66 for(int i=0;i<=len;i++) cop[i].x/=len; 67 } 68 void Square(cpx *cop,int len) 69 { 70 register int i; 71 Trans(cop,len,1); 72 for(i=0;i<=len;i++) 73 cop[i]=cop[i]*cop[i]; 74 Trans(cop,len,-1); 75 } 76 double Ask(double *f,double x,int len) 77 { 78 double ret=0,var=1; 79 for(int i=0;i<=len;i++) 80 ret+=f[i]*var,var*=x; 81 return ret; 82 } 83 double Iterate(double x) 84 { 85 for(int i=1;i<=30;i++) 86 { 87 double y=Ask(f,x,n); 88 if(fabs(y)<=eps) return x; 89 x-=y/Ask(d,x,nm),x=max(x,ll),x=min(x,rr); 90 } 91 return inf; 92 } 93 int main() 94 { 95 register int i; 96 scanf("%d%d%d%lf%lf",&la,&lb,&lc,&ll,&rr); 97 for(i=0;i<=la;i++) Read(a[i].x); 98 for(i=0;i<=lb;i++) Read(b[i].x); 99 for(i=0;i<=lc;i++) Read(c[i].x); 100 Prework(),Square(a,n),Square(b,n),Square(c,n); 101 for(i=0;i<=n;i++) c[i]=c[i]-a[i]-b[i]; nm=n-1; 102 for(i=0;i<=n;i++) f[i]=c[i].x,d[i]=(i+1)*c[i+1].x; 103 double ans=Iterate((ll+rr)/2); 104 ans>=inf?printf("Inconsistent!"):printf("%.8f",ans); 105 return 0; 106 }