A.Rational Approximation回到顶部
题意
给一个m+n(10)位多项式f=f0+f1x+f2x^2...,求一个最多n-1位多项式q和最多m-1位多项式p,使得q*f-p得到的最低项的指数为n+m-1,系数为1
题解
大概就是设一个多项式q,然后q*f,高斯消元解出q,最后计算得到p,复杂度O(10^3)
代码
1 #include<stdio.h> 2 #include<algorithm> 3 #include<iostream> 4 #include<queue> 5 #include<string.h> 6 using namespace std; 7 8 #define LL long long 9 class rat 10 { 11 public: 12 LL num,den; 13 rat() { 14 num=0;den=1; 15 }; 16 rat(LL n,LL d){ 17 bool pos=true; 18 if(d<0){ 19 d=-d;pos=!pos; 20 } 21 if(n<0){ 22 n=-n;pos=!pos; 23 } 24 LL fac=__gcd(n,d); 25 num=n/fac; 26 den=d/fac; 27 if(!pos)num=-num; 28 }; 29 rat add(rat rhs){ 30 LL newNum=num*rhs.den+den*rhs.num; 31 LL newDen=den*rhs.den; 32 rat ans(newNum,newDen); 33 return ans; 34 }; 35 rat sub(rat rhs){ 36 LL newNum=num*rhs.den-den*rhs.num; 37 LL newDen=den*rhs.den; 38 rat ans(newNum,newDen); 39 return ans; 40 }; 41 rat mul(rat rhs){ 42 LL newNum=num*rhs.num; 43 LL newDen=den*rhs.den; 44 rat ans(newNum,newDen); 45 return ans; 46 }; 47 rat div(rat rhs){ 48 LL newNum=num*rhs.den; 49 LL newDen=den*rhs.num; 50 rat ans(newNum,newDen); 51 return ans; 52 }; 53 void print(){ 54 if(num==0)cout<<'0'; 55 else if(den==1)cout<<num; 56 else cout<<num<<'/'<<den; 57 } 58 }; 59 60 rat operator+(rat a,rat b){return a.add(b);}; 61 rat operator-(rat a,rat b){return a.sub(b);}; 62 rat operator*(rat a,rat b){return a.mul(b);}; 63 rat operator/(rat a,rat b){return a.div(b);}; 64 ostream& operator<<(ostream &out, rat a){a.print();return out;}; 65 66 const int N=15; 67 68 rat q[N],p[N],f[N],a[N][N],b[N]; 69 70 int n,m; 71 void printarray()//debug 72 { 73 for(int i=0; i<n; i++) { 74 for(int j=0; j<n; j++) 75 cout << a[i][j] << ' '; 76 cout<<b[i]; 77 cout << endl; 78 } 79 cout<<endl; 80 } 81 void Gauss(){ 82 for(int i=0;i<n-1;i++){ 83 int j=i; 84 while(a[i][j].num==0)j++; 85 if(j!=i){ 86 rat temp=b[j];b[j]=b[i];b[i]=temp; 87 for(int k=0;k<n;k++){ 88 temp=a[j][k];a[j][k]=a[i][k];a[i][k]=temp; 89 } 90 } 91 for(int k=i+1;k<n;k++){ 92 rat fact=a[k][i]/a[i][i]; 93 for(int l=i;l<n;l++){ 94 a[k][l]=a[k][l]-a[i][l]*fact; 95 } 96 b[k]=b[k]-b[i]*fact; 97 } 98 //printarray(); 99 } 100 q[n-1]=b[n-1]/a[n-1][n-1]; 101 for(int i=n-2;i>=0;i--){ 102 q[i]=b[i]; 103 for(int j=i+1;j<n;j++)q[i]=q[i]-a[i][j]*q[j]; 104 q[i]=q[i]/a[i][i]; 105 } 106 } 107 void pr(rat *s,int r){ 108 int out=0; 109 for(int i=0;i<r;i++){ 110 if(s[i].num==0)continue; 111 if(out==0)out=1; 112 else cout<<" "; 113 cout<<"("<<s[i]<<","<<i<<")"; 114 } 115 if(out==0)cout<<"(0,0)"; 116 cout<<'\n'; 117 } 118 int main(){ 119 int ca=0; 120 while(cin>>m>>n,m||n){ 121 if(ca++)cout<<'\n'; 122 for(int i=0;i<m+n;i++)cin>>f[i].num,f[i].den=1; 123 //n*n 124 //a x^n+m-1 x^n+m-2 .... x^0 125 //b x^n+m-1 x^n+m-2 .... x^0 126 for(int i=0;i<n;i++){ 127 for(int j=0;j<n;j++){ 128 if(m+j-i>=0){ 129 //printf("i=%d n-j-1=%d f[%d]=%lld\n",i,n-j-1,m+j-i,f[m+j-i]); 130 a[i][n-j-1]=f[m+j-i]; 131 } 132 else{ 133 //printf("i=%d n-j-1=%d 0 1\n",i,n-j-1); 134 a[i][n-j-1].num=0; 135 a[i][n-j-1].den=1; 136 } 137 } 138 b[i].num=(i==0)?1:0; 139 b[i].den=1; 140 //printf("i=%d num=%lld\n",i,b[i]); 141 } 142 Gauss(); 143 //for(int i=0;i<n;i++) 144 // q[i].print(),printf(" "); 145 for(int i=0;i<m;i++)p[i].num=0,p[i].den=1; 146 for(int i=0;i<n+m;i++) 147 for(int j=0;j<n;j++) 148 if(i+j<m) 149 p[i+j]=p[i+j]+f[i]*q[j]; 150 pr(p,m);pr(q,n); 151 } 152 return 0; 153 }