zoj 2369 Two Cylinders
链接：http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=2369
题意：已知两个无限长的圆柱的半径，它们轴线相交且垂直，求相交部分的体积。
思路：以相交点为坐标原点，建立三维坐标系，列出 两个圆柱的公式，把三维转化为一维积分公式
V1 : x2 + y2 = r12 ==> y = sqrt ( r12 - x2 );
V2 : x2 + z2 = r22 ==> z = sqrt ( r22 - x2 );
V = ∫∫∫ dv = ∫∫∫ dxdydz = ∫(∫∫ dydz)dx 积分区间：( 0 <= x <= r1 , 0 <= x <= r2 )
所以： V = ∫0min(r1, r2) sqrt(( r12-x2) * (r22-x2));
计算这个积分，我用了好几种方法，但是只有 变步长simpon 积分达到精度要求，我也不清楚什么原因。下面是几种求积分的算法：
测试数据：
10 1 1 2 1 3 2 5 2 8 3 8 4 99 2 99 56 84 31 91 17
正确答案：
5.3333 12.1603 70.9361 123.0975 444.2909 778.2605 2488.0144 1869197.4256 498415.2860 164517.4621
变步长 simpon 积分：（得到正确结果）
1 #include <iostream> 2 #include <stdio.h> 3 #include <algorithm> 4 #include <queue> 5 #include <map> 6 #include <set> 7 #include <vector> 8 #include <cstring> 9 #include <math.h> 10 using namespace std; 11 12 const double eps = 1e-4; 13 double r1, r2; 14 15 double f(double x) 16 { 17 return 8. * sqrt((r1*r1 - x*x) * (r2*r2 - x*x)); 18 } 19 20 double simpson(double a, double b) 21 { 22 int n, k; 23 double h, t1, t2, s1, s2, ep, p, x; 24 n = 1, h = b - a; 25 26 t1 = h * ( f(a) + f(b) ) / 2.; //计算T1 = (b-a)/2*[f(a) + f(b)] 27 s1 = t1; //T1代替S1 28 ep = eps + 1.; 29 while(ep >= eps) 30 { 31 p = 0; 32 for(k = 0; k <= n-1; k++) 33 { 34 x = a + (k + 0.5) * h; 35 p += f(x); 36 } 37 t2 = (t1 + h * p) / 2.; //计算T2n = Tn/2 + h/2 * [(累加)f(xk + 0.5)]; 38 s2 = (4. * t2 - t1) / 3.; //计算Sn = (4*Tn - Tn) / 3; 39 ep = fabs(s2 - s1); //计算精度 40 t1 = t2, s1 = s2, n += n, h /= 2.; 41 } 42 return s2; 43 } 44 45 int main() 46 { 47 int CASE; 48 freopen("zoj2369.txt", "r", stdin); 49 cin >> CASE; 50 while(CASE--) 51 { 52 scanf("%lf%lf", &r1, &r2); 53 if(r2 < r1) swap(r1, r2); 54 double t = simpson(0., r1); 55 printf("%.4lf ", t); 56 57 } 58 return 0; 59 } 60 61 62 /* 结果： 63 5.3333 64 12.1603 65 70.9361 66 123.0974 67 444.2908 68 778.2604 69 2488.0144 70 1869197.4256 71 498415.2860 72 164517.4621 73 */
勒让德-高斯求积（精度不够）：
1 #include <iostream> 2 #include <stdio.h> 3 #include <algorithm> 4 #include <queue> 5 #include <map> 6 #include <set> 7 #include <vector> 8 #include <cstring> 9 #include <math.h> 10 using namespace std; 11 12 const double eps = 1e-4; 13 double r1, r2; 14 15 double f(double x) 16 { 17 return 8. * sqrt((r1*r1 - x*x) * (r2*r2 - x*x)); 18 } 19 20 double lrgs(double a, double b) 21 { 22 int m, i, j; 23 double s, p, ep, aa, bb, w, x, g, h; 24 static double t[5] = {-0.9061798459, -0.5384693101, 0.0, 0.5384693101, 0.9061798459}; 25 static double c[5] = {0.2369268851, 0.4786286705, 0.5688888889, 0.4786286705, 0.2369268851}; 26 m = 1; 27 h = b - a; s = fabs(0.001 * h); 28 p = 1.0e+35, ep = eps + 1.0; 29 while(ep >= eps && fabs(h) > s) 30 { 31 g = 0.; 32 for(i = 1; i <= m; i++) 33 { 34 aa = a + (i-1.) * h, bb = a + i * h; 35 w = 0.; 36 for(j = 0; j < 5; j++) 37 { 38 x = ((bb - aa) * t[j] + (bb + aa)) / 2.; 39 w += f(x) * c[j]; 40 } 41 g += w; 42 } 43 g = g * h / 2.; 44 ep = fabs(g-p) / (1. + fabs(g)); 45 p = g, m += 1, h = (b-a) / m; 46 } 47 return g; 48 } 49 50 int main() 51 { 52 int t; 53 freopen("zoj2369.txt", "r", stdin); 54 scanf("%d", &t); 55 while(t--) 56 { 57 scanf("%lf%lf", &r1, &r2); 58 if(r1 > r2) swap(r1, r2); 59 printf("%.4lf ", lrgs(0, r1)); 60 } 61 return 0; 62 } 63 64 /*结果 65 5.3333 66 12.1619 67 70.9440 68 123.1138 69 444.3504 70 778.3592 71 2488.3679 72 1869425.2091 73 498482.1853 74 164540.5195 75 */
龙贝格求积法（精度不够）：
1 #include <iostream> 2 #include <stdio.h> 3 #include <algorithm> 4 #include <vector> 5 #include <cstring> 6 #include <math.h> 7 using namespace std; 8 9 const double eps = 1e-4; 10 double r1, r2; 11 12 double f(double x) 13 { 14 return 8. * sqrt((r1*r1 - x*x) * (r2*r2 - x*x)); 15 } 16 17 double romb(double a, double b) 18 { 19 int n, m, i, k; 20 double y[10], h, ep, p, x, s, q; 21 h = b - a; 22 y[0] = h * ((f(a) + f(b))) / 2.; 23 m = 1, n = 1, ep = eps + 1.0; 24 while(ep >= eps ) 25 { 26 p = 0.0; 27 for(i = 0; i < n; i++) 28 { 29 x = a + (i+0.5) * h; 30 p += f(x); 31 } 32 p = (y[0] + h * p) /2.; 33 s = 1.0; 34 for(k = 1; k <= m; k++) 35 { 36 s = 4. * s; 37 q = (s*p - y[k-1]) / (s - 1.); 38 y[k - 1] = p, p = q; 39 } 40 ep = fabs(q - y[m-1]); 41 m = m + 1, y[m - 1] = q, n += n, h /= 2.; 42 } 43 return q; 44 } 45 46 int main() 47 { 48 int t; 49 freopen("zoj2369.txt", "r", stdin); 50 freopen("zzin1.txt", "w", stdout); 51 scanf("%d", &t); 52 while(t--) 53 { 54 scanf("%lf%lf", &r1, &r2); 55 if(r1 > r2) swap(r1, r2); 56 printf("%.4lf ", romb(0, r1)); 57 } 58 return 0; 59 } 60 61 /*结果 62 5.3333 63 12.1529 64 70.8978 65 123.0697 66 444.1897 67 778.0925 68 2487.8023 69 1869191.3654 70 498410.2627 71 164515.7323 72 */
也不知最后两个代码那里写搓了，真心无力。
下面给出标程：
1 #include <cmath> 2 #include <cstdio> 3 #include <algorithm> 4 5 using namespace std; 6 7 const long double H = 5e-5; 8 9 struct F { 10 long double rr1, rr2; 11 F(long double r1, long double r2) : rr1(r1 * r1), rr2(r2 * r2) { 12 } 13 long double operator()(long double x) const { 14 x = x * x; 15 return sqrt((rr2 - x) * (rr1 - x)); 16 } 17 }; 18 19 int main() { 20 int re; 21 22 scanf("%d", &re); 23 for (int ri = 1; ri <= re; ++ri) { 24 long double r1, r2; 25 scanf("%Lf%Lf", &r1, &r2); 26 if (r1 > r2) { 27 swap(r1, r2); 28 } 29 int n = int(r1 / H); 30 // int n = 1000000; 31 long double h = r1 / n / 2; 32 F f(r1, r2); 33 long double s = 0; 34 s += f(0) + f(r1); 35 for (int i = 1; i < 2 * n; ++i) { 36 if (i & 1) { 37 s += 4 * f(i * h); 38 } else { 39 s += 2 * f(i * h); 40 } 41 } 42 s /= 6 * n; 43 s *= 8 * r1; 44 printf("%.6Lf ", s); 45 } 46 47 return 0; 48 }
网上看的大牛的代码，0ms秒过，我的用了290ms， 标程用了310ms
1 #include <cstdio> 2 #include <cmath> 3 #include <algorithm> 4 using namespace std; 5 6 #define db double 7 #define rt return 8 #define cs const 9 10 cs db eps = 1e-9; 11 inline int sig(db x){rt (x>eps)-(x<-eps);} 12 13 db r1, r2; 14 15 inline db f(db x) { 16 rt 4. * sqrt(r1*r1-x*x) * sqrt(r2*r2-x*x); 17 } 18 19 inline db simpson(db a, db b) { 20 rt (b-a)*(f(a) + 4.*f((a+b)/2.) + f(b)) / 6.; 21 } 22 23 inline db calc(db a, db b, int depth) { 24 db total = simpson(a, b), mid = (a+b) / 2.; 25 db tmp = simpson(a, mid) + simpson(mid, b); 26 if(sig(total-tmp) == 0 && depth > 3) rt total; 27 rt calc(a, mid, depth + 1) + calc(mid, b, depth + 1); 28 } 29 30 int main() { 31 // freopen("std.in", "r", stdin); 32 int T; 33 freopen("zoj2369.txt", "r", stdin); 34 scanf("%d", &T); 35 while(T--) { 36 scanf("%lf%lf", &r1, &r2); 37 if(sig(r1-r2) > 0) swap(r1, r2); 38 printf("%.4lf ", 2. * calc(0., r1, 0)); 39 } 40 rt 0; 41 42 }