2016北京集训测试赛(十六)Problem C: ball

2016北京集训测试赛(十六)Problem C: ball

Solution

这是一道好题.

考虑球体的体积是怎么计算的: 我们令\(f_k(r)\)表示\(x\)维单位球的体积, 则

\[f_k(1) = \int_{-1}^1 f_{k - 1}(\sqrt{1 - r^2}) dr
\]

然而\(f_{k - 1}(\sqrt{1 - x^2})\)并不容易处理, 我们又注意到\(k\)维球体的体积可以表示为\(a \pi r^k\), 因此\(f_k(\sqrt{1 - r^2}) = f_k(1) \times (1 - r)^{\frac k 2}\)

因此递归式变成了这个样子:

\[f_k(1) = \int_{-1}^1 f_{k - 1}(1) \times (1 - r)^{\frac{k - 1}2} dr = f_{k - 1} \int_{-1}^1 (1 - r)^{\frac{k - 1}2} dr
\]

Simpson积分出每个\(f_k(1)\), 递推即可.

哦对了, 差点忘了最后要进行的线性变换. 直接在原来体积的基础上乘上该矩阵行列式的值即可.

#include <cstdio>
#include <algorithm>
#include <cmath>
#define swap std::swap const int D = 19;
const double EPS = 1e-12;
struct determinant
{
double a[D][D];
inline double elemination(int n)
{
double ans = 1;
for(int i = 0; i < n; ++ i)
{
int p = i; for(; p < n && a[p][i] == 0; ++ p);
if(p != i) for(int j = 0; j < n; ++ j) swap(a[i][j], a[p][j]);
for(int j = i + 1; j < n; ++ j) if(a[j][i] != 0)
{
double tmp = a[j][i] / a[i][i];
for(int k = 0; k < n; ++ k) a[j][k] -= tmp * a[i][k];
}
}
for(int i = 0; i < n; ++ i) ans *= a[i][i];
return fabs(ans);
}
}det;
inline double calculate(double a, double x)
{
double res = 1; for(int i = 0; i < x - 1; ++ i) res *= 1 - a * a; return sqrt(res);
}
inline double integrate(double L, double R, double x)
{
double a = calculate(L, x), b = calculate(R, x), c = calculate((L + R) / 2, x);
if(fabs(c * (R - L) - (a + b) * (R - L) / 2) > EPS) return integrate(L, (L + R) / 2, x) + integrate((L + R) / 2, R, x);
else return c * (R - L);
}
inline double work(int d)
{
if(d == 1) return 2;
double lst = work(d - 1);
return lst * integrate(-1, 1, d);
}
int main()
{ #ifndef ONLINE_JUDGE freopen("ball.in", "r", stdin);
freopen("ball.out", "w", stdout); #endif int d; scanf("%d", &d);
for(int i = 0; i < d; ++ i) for(int j = 0; j < d; ++ j) scanf("%lf", &det.a[i][j]);
printf("%.10lf\n", work(d) * det.elemination(d));
}
上一篇:2016北京集训测试赛(十)Problem A: azelso


下一篇:uva----11729 Commando war (突击战争)