LA 6135 - Environment Protection simpson积分

求深度 d 使得从深度0到深度d 的middle layer的面积是 A。

一开始读错题意。 所以代码中的 A 是深度 ，ans 才是面积A。 变量名不一致。

#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <stack>
#include <cstring>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <assert.h>
#include <queue>
#define REP(i,n) for(int i=0;i<n;i++)
#define TR(i,x) for(typeof(x.begin()) i=x.begin();i!=x.end();i++)
#define ALLL(x) x.begin(),x.end()
#define SORT(x) sort(ALLL(x))
#define CLEAR(x) memset(x,0,sizeof(x))
#define FILLL(x,c) memset(x,c,sizeof(x))
using namespace std;
const double eps = 1e-9;
#define LL long long 
#define pb push_back
const int maxn  = 11000;

int W,D,k;
double A,ans;
double a[15],b[15],c[15],d[15];
int flag ;
double F(double x){  
      if(flag==1){
      	   double up = a[0];
      	   double tmp = x;
      	   for(int i=1;i<=k;i++){
      	        	up+= a[i]*tmp;
      	        	tmp *= x;
      	   }
      	   double down = b[0];
      	   tmp = x;
      	   for(int i=1;i<=k;i++){
      	   	    down += b[i]*tmp;
      	   	    tmp *= x;
      	   }
      	   //cout << up<< "  uppppp "<<down <<"  "<<x<<endl;
      	   double ans = up/down ;
      	   
      	   return  max(ans,-(A+0.0));
      } else if(flag ==2){
      	   double up = c[0];
      	   double tmp = x;
      	   for(int i=1;i<=k;i++){
      	        	up+= c[i]*tmp;
      	        	tmp *= x;
      	   }
      	   double down = d[0];
      	   tmp = x;
      	   for(int i=1;i<=k;i++){
      	   	    down += d[i]*tmp;
      	   	    tmp *= x;
      	   }
      	  // cout << up << "  "<<down<<endl;
      	   double ans = up/down ;
      	   return  max(ans,-(A+0.0));
      } 
      return 0;
}  
  
// 三点simpson法。这里要求F是一个全局函数    
double simpson(double a,double b){  
    double c =  a+(b-a)/2;  
    return (F(a) + 4*F(c) + F(b))*(b-a)/6;  
}  
// 自适应Simpson公式（递归过程）。已知整个区间[a,b]上的三点simpson值A   
double asr(double a , double b ,double eps ,double A){  
//	cout << a << " "<< b<< " "<<endl;
    double c = a+ (b-a)/2;  
    double L = simpson(a,c) ,R = simpson(c,b);  
   // cout << A << "  "<< L<< "  "<<R<<endl;
    if(fabs(A-L-R)<=15*eps) return L + R +(A-L-R)/15;  
    return asr(a,c,eps/2,L) + asr(c,b,eps/2,R);  
}  
// 自适应Simpson公式（主过程）    
double asr(double a, double b, double eps) {    
  return asr(a, b, eps, simpson(a, b));    
} 

double  get_ans(){
	flag =1;
	double ans = asr(0,W+0.0,1e-6);
	//cout << ans <<endl;
    flag =2 ;
    double ans2 = asr(0,W+0.0,1e-8);
    //cout << ans2 <<endl;
	ans -= ans2;
	return ans;
}
void solve(){
	 double left = 0;
	 double right = D+0.0;
	 int t = 100;
	 while(right-left>1e-9){
	     double  mid = (left + right)/2;
	     A = mid;
	     double t = get_ans();
	     if(t<ans){
	     	  left = mid;
	     }else{
	     	   right = mid;
	     }
	 }
	 printf("%.5lf\n",left);
	
}
int main(){
     while(~scanf("%d%d",&W,&D)){
     	scanf("%lf",&ans);
     	scanf("%d",&k);
     	for(int i=0;i<=k;i++)scanf("%lf",&a[i]);
     	for(int i=0;i<=k;i++)scanf("%lf",&b[i]);
     	for(int i=0;i<=k;i++)scanf("%lf",&c[i]);
     	for(int i=0;i<=k;i++)scanf("%lf",&d[i]);
     	solve();
     }
    return 0;
}

求F 的函数更改。。 积分不再是求两次

时间快了不少

#include <vector>  
#include <list>  
#include <map>  
#include <set>  
#include <deque>  
#include <stack>  
#include <cstring>  
#include <bitset>  
#include <algorithm>  
#include <functional>  
#include <numeric>  
#include <utility>  
#include <sstream>  
#include <iostream>  
#include <iomanip>  
#include <cstdio>  
#include <cmath>  
#include <cstdlib>  
#include <ctime>  
#include <assert.h>  
#include <queue>  
#define REP(i,n) for(int i=0;i<n;i++)  
#define TR(i,x) for(typeof(x.begin()) i=x.begin();i!=x.end();i++)  
#define ALLL(x) x.begin(),x.end()  
#define SORT(x) sort(ALLL(x))  
#define CLEAR(x) memset(x,0,sizeof(x))  
#define FILLL(x,c) memset(x,c,sizeof(x))  
using namespace std;  
const double eps = 1e-9;  
#define LL long long   
#define pb push_back  
const int maxn  = 11000;  
  
int W,D,k;  
double A,ans;  
double a[15],b[15],c[15],d[15];  
int flag ;  
double F(double x){    
      
           double up = a[0];  
           double tmp = x;  
           for(int i=1;i<=k;i++){  
                    up+= a[i]*tmp;  
                    tmp *= x;  
           }  
           double down = b[0];  
           tmp = x;  
           for(int i=1;i<=k;i++){  
                down += b[i]*tmp;  
                tmp *= x;  
           }  
           //cout << up<< "  uppppp "<<down <<"  "<<x<<endl;  
           double ans = up/down ;  
             
           ans = max(ans,-(A+0.0));  
       
           double up2 = c[0];  
           double tmp2 = x;  
           for(int i=1;i<=k;i++){  
                    up2+= c[i]*tmp2;  
                    tmp2 *= x;  
           }  
           double down2 = d[0];  
           tmp2 = x;  
           for(int i=1;i<=k;i++){  
                down2 += d[i]*tmp2;  
                tmp2 *= x;  
           }  
          // cout << up << "  "<<down<<endl;  
           double ans3 = up2/down2 ;  
           ans3 =  max(ans3,-(A+0.0));  
         
      return ans- ans3;  
}    
    
// 三点simpson法。这里要求F是一个全局函数      
double simpson(double a,double b){    
    double c =  a+(b-a)/2;    
    return (F(a) + 4*F(c) + F(b))*(b-a)/6;    
}    
// 自适应Simpson公式（递归过程）。已知整个区间[a,b]上的三点simpson值A     
double asr(double a , double b ,double eps ,double A){    
//  cout << a << " "<< b<< " "<<endl;  
    double c = a+ (b-a)/2;    
    double L = simpson(a,c) ,R = simpson(c,b);    
   // cout << A << "  "<< L<< "  "<<R<<endl;  
    if(fabs(A-L-R)<=15*eps) return L + R +(A-L-R)/15;    
    return asr(a,c,eps/2,L) + asr(c,b,eps/2,R);    
}    
// 自适应Simpson公式（主过程）      
double asr(double a, double b, double eps) {      
  return asr(a, b, eps, simpson(a, b));      
}   
  
double  get_ans(){  
     
    double ans = asr(0,W+0.0,1e-6);  
    //cout << ans <<endl;  
  
    return ans;  
}  
void solve(){  
     double left = 0;  
     double right = D+0.0;  
     int t = 100;  
     while(right-left>1e-9){  
         double  mid = (left + right)/2;  
         A = mid;  
         double t = get_ans();  
         if(t<ans){  
              left = mid;  
         }else{  
               right = mid;  
         }  
     }  
     printf("%.5lf\n",left);  
      
}  
int main(){  
     while(~scanf("%d%d",&W,&D)){  
        scanf("%lf",&ans);  
        scanf("%d",&k);  
        for(int i=0;i<=k;i++)scanf("%lf",&a[i]);  
        for(int i=0;i<=k;i++)scanf("%lf",&b[i]);  
        for(int i=0;i<=k;i++)scanf("%lf",&c[i]);  
        for(int i=0;i<=k;i++)scanf("%lf",&d[i]);  
        solve();  
     }  
    return 0;  
}