Simpson's rule - Wikipedia, the free encyclopedia
利用这个公式，用二分的方法来计算积分。
1 #include <iostream> 2 #include <algorithm> 3 #include <cstdio> 4 #include <cstring> 5 #include <cmath> 6 7 using namespace std; 8 9 const double EPS = 1e-8; 10 double A, B, C, P, Q; 11 12 template<class T> T sqr(T x) { return x * x;} 13 inline double cal(double x) { return A * sqr(x) + (B - P) * x + C - Q;} 14 inline double sps(double l, double r) { return (cal(l) + cal(r) + 4 * cal((l + r) / 2)) / 6 * (r - l);} 15 16 double work(double l, double r) { 17 //cout << l << ' ' << r << endl; 18 double ans = sps(l, r), m = (l + r) / 2; 19 if (fabs(ans - sps(l, m) - sps(m, r)) < EPS) return ans; 20 else return work(l, m) + work(m, r); 21 } 22 23 24 int main() { 25 int T; 26 double l, r; 27 double x[3], y[3]; 28 cin >> T; 29 while (T--) { 30 for (int i = 0; i < 3; i++) cin >> x[i] >> y[i]; 31 double p[2], q[2], d[2]; 32 for (int i = 0; i < 2; i++) p[i] = sqr(x[i]) - sqr(x[i + 1]), q[i] = x[i] - x[i + 1], d[i] = y[i] - y[i + 1]; 33 A = (q[1] * d[0] - q[0] * d[1]) / (p[0] * q[1] - p[1] * q[0]); 34 B = (p[1] * d[0] - p[0] * d[1]) / (p[1] * q[0] - p[0] * q[1]); 35 C = y[0] - B * x[0] - A * sqr(x[0]); 36 //cout << A << ' ' << B << ' ' << C << endl; 37 P = (y[1] - y[2]) / (x[1] - x[2]); 38 Q = y[1] - P * x[1]; 39 //cout << P << ' ' << Q << endl; 40 printf("%.2f ", work(x[1], x[2])); 41 } 42 return 0; 43 }
1 #include <iostream> 2 #include <algorithm> 3 #include <cstdio> 4 #include <cstring> 5 #include <cmath> 6 7 using namespace std; 8 9 const double EPS = 1e-8; 10 double A, B; 11 12 template<class T> T sqr(T x) { return x * x;} 13 inline double cal(double x) { return 2 * B * sqrt(1 - sqr(x) / sqr(A));} 14 inline double sps(double l, double r) { return (cal(l) + cal(r) + 4 * cal((l + r) / 2)) / 6 * (r - l);} 15 16 double work(double l, double r) { 17 //cout << l << ' ' << r << endl; 18 double ans = sps(l, r), m = (l + r) / 2; 19 if (fabs(ans - sps(l, m) - sps(m, r)) < EPS) return ans; 20 else return work(l, m) + work(m, r); 21 } 22 23 int main() { 24 int T; 25 double l, r; 26 cin >> T; 27 while (T-- && cin >> A >> B >> l >> r) printf("%.3f ", work(l, r)); 28 return 0; 29 }
之后还有题会继续更新。
UPD：
就是因为见过这题，所以才学这个公式的。1y~
1 #include <cstdio> 2 #include <cstring> 3 #include <iostream> 4 #include <algorithm> 5 #include <cmath> 6 7 using namespace std; 8 9 double coe[4][11]; 10 const double EPS = 1e-8; 11 12 int k; 13 double cal(double x, double *c) { 14 double ret = c[0]; 15 for (int i = 1; i <= k; i++) ret *= x, ret += c[i]; 16 return ret; 17 } 18 19 inline double cal(double x, double *p, double *q) { return cal(x, p) / cal(x, q);} 20 inline double cal(double x, double y, double *p, double *q) { return max(cal(x, p, q) - y, 0.0);} 21 inline double simpson(double y, double l, double r, double *p, double *q) { return (cal(l, y, p, q) + cal(r, y, p, q) + 4 * cal((l + r) / 2, y, p, q)) * (r - l) / 6;} 22 23 inline double getpart(double y, double l, double r, double *p, double *q) { 24 double sum = simpson(y, l, r, p, q); 25 //cout << l << ' ' << r << ' ' << sum << endl; 26 if (fabs(sum - simpson(y, l, (l + r) / 2, p, q) - simpson(y, (l + r) / 2, r, p, q)) < EPS) return sum; 27 return getpart(y, l, (l + r) / 2, p, q) + getpart(y, (l + r) / 2, r, p, q); 28 } 29 30 inline double getarea(double y, double l, double r, double *p, double *q) { 31 double ret = 0, d = (r - l) / 100; 32 for (int i = 0; i < 100; i++) { 33 ret += getpart(y, l + d * i, l + d * (i + 1), p, q); 34 } 35 return ret; 36 } 37 38 double dc2(double l, double r, double a, double w) { 39 double m; 40 while (r - l > EPS) { 41 m = (l + r) / 2.0; 42 //cout << m << ' ' << getarea(m, 0, w, coe[0], coe[1]) - getarea(m, 0, w, coe[2], coe[3]) << endl; 43 if (getarea(m, 0, w, coe[0], coe[1]) - getarea(m, 0, w, coe[2], coe[3]) > a) l = m; 44 else r = m; 45 } 46 return l; 47 } 48 49 int main() { 50 //freopen("in", "r", stdin); 51 //freopen("out", "w", stdout); 52 double w, d, a; 53 while (cin >> w >> d >> a >> k) { 54 for (int i = 0; i < 4; i++) for (int j = 0; j <= k; j++) cin >> coe[i][j]; 55 for (int i = 0; i < 4; i++) reverse(coe[i], coe[i] + k + 1); 56 //cout << getarea(-5.51389, 0, w, coe[0], coe[1]) - getarea(-5.51389, 0, w, coe[2], coe[3]) << endl; 57 //cout << cal(3, coe[0], coe[1]) << endl; 58 printf("%.5f ", -dc2(-d, 0, a, w)); 59 } 60 return 0; 61 }
——written by Lyon