在一个正方形内随机选取一点, 并将此点与四个顶点连上直线, 从而将正方形分割为四个三角形. 求四个三角形之中所有内角均不超过 120° 的概率, 答案准确至小数点后 3 位.
@tvb 在她的方案中 讲了的数学方法需要用到各种三角公式来计算面积, 我数学功底不行, 这种方法对我来说太麻烦(虽然算面积这种想法我也想到了).
暴力, 我想大家都会首先想到, 而且程序好写, 需要注意的一点就是你起始点不能在正方形边上, 然后 x 和 y 坐标 都以一个非常小的 步长 增量, 比如 x += 1e-5, y += 1e-5.
补一句: 这里的暴力并不是真的暴力, 因为你无法穷举所有的点, 这里只是尽可能的穷举“所有”, 有一点小误差.
这样一来, 程序的思路就非常明确了, 只需要判断一下是否是有效点, 然后统计一下除于总个数就是概率.
判断逻辑也较为简单, 根据 (x, y) 连接四个点, 然后得到4个三角形, 依次对每个三角形进行判断内角, 这里我们需要用到一个三角公式.
已知 三边长, 可以通过 公式: c^2=a^2+b^2−2abcos(C) 来推导三个角的 cos 值, 然后你可以 反余玄来获取角度, 但也可以直接 判断 是否大于 -0.5. 因为 cos(120度)=-0.5, 余玄值若是大于 -0.5 则角度小于 120度.
假设给定两个点的坐标, 我们算出边长:
1 2 3 | function getSide(x1, y1, x2, y2) { return Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) } |
function getSide(x1, y1, x2, y2) { return Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) }
所以判断逻辑为:
1 2 3 4 5 6 7 8 9 | function check(x1, y1, x2, y2, x3, y3) { var a = getSide(x1, y1, x2, y2); var b = getSide(x1, y1, x3, y3); var c = getSide(x2, y2, x3, y3); var cosc = (a * a + b * b - c * c) / (2 * a * b); var cosa = (b * b + c * c - a * a) / (2 * b * c); var cosb = (c * c + a * a - b * b) / (2 * c * a); return (cosc >= -0.5) && (cosb >= -0.5) && (cosa >= -0.5); } |
function check(x1, y1, x2, y2, x3, y3) { var a = getSide(x1, y1, x2, y2); var b = getSide(x1, y1, x3, y3); var c = getSide(x2, y2, x3, y3); var cosc = (a * a + b * b - c * c) / (2 * a * b); var cosa = (b * b + c * c - a * a) / (2 * b * c); var cosb = (c * c + a * a - b * b) / (2 * c * a); return (cosc >= -0.5) && (cosb >= -0.5) && (cosa >= -0.5); }
如果你是用 Math.acos 算出角度, 可以变成:
1 2 3 4 5 6 7 8 9 10 11 | function check(x1, y1, x2, y2, x3, y3) { var a = getSide(x1, y1, x2, y2); var b = getSide(x1, y1, x3, y3); var c = getSide(x2, y2, x3, y3); var cosc = (a * a + b * b - c * c) / (2 * a * b); var cosa = (b * b + c * c - a * a) / (2 * b * c); angle_c = Math.acos(cosc); angle_a = Math.acos(cosa); angle_b = 180 - angle_c - angle_a; return (angle_c <= 120) && (angle_b <= 120) && (angle_a <= 120); } |
function check(x1, y1, x2, y2, x3, y3) { var a = getSide(x1, y1, x2, y2); var b = getSide(x1, y1, x3, y3); var c = getSide(x2, y2, x3, y3); var cosc = (a * a + b * b - c * c) / (2 * a * b); var cosa = (b * b + c * c - a * a) / (2 * b * c); angle_c = Math.acos(cosc); angle_a = Math.acos(cosa); angle_b = 180 - angle_c - angle_a; return (angle_c <= 120) && (angle_b <= 120) && (angle_a <= 120); }
接下来撒点就是水到渠成的事了:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | var total = 0; // 所有点数 var valid = 0; // 符合条件的点 for (let x = 1e-7; x <= 0.99999; x += 1e-4) { for (let y = 1e-7; y <= 0.99999; y += 1e-4) { total ++; if ( check(x, y, 0, 0, 0, 1) && check(x, y, 0, 1, 1, 1) && check(x, y, 1, 1, 1, 0) && check(x, y, 1, 0, 0, 0) ) { valid ++; } } } console.log(valid); console.log(total); console.log(Math.round(valid / total * 1000) / 1000); // 保留三位精度 |
var total = 0; // 所有点数 var valid = 0; // 符合条件的点 for (let x = 1e-7; x <= 0.99999; x += 1e-4) { for (let y = 1e-7; y <= 0.99999; y += 1e-4) { total ++; if ( check(x, y, 0, 0, 0, 1) && check(x, y, 0, 1, 1, 1) && check(x, y, 1, 1, 1, 0) && check(x, y, 1, 0, 0, 0) ) { valid ++; } } } console.log(valid); console.log(total); console.log(Math.round(valid / total * 1000) / 1000); // 保留三位精度
增量试着从 1e-3 到 1e-6, 因为小数点需要3位, 所以 1e-3 之后精度就没有再变化了(天啊撸, 我提交答案的时候以为需要的三位数是(乘以100后的) 21.255% 所以花了不少时间做了无用工, 而且提交上去答案还算是错的了)
以上方法, 不管你跑多少次, 输出的结果一定是一样的, 计算机学上我们称这种算法为 确定型 (Deterministic)
这种方法的复杂度是 O(xy), 当步长很小的时候需要的时候很不环保, 我们可以试着用蒙特卡罗防真的方法: 当我们随机往正方形内撒点的时候, 撒足够多的点, 那么符合条件的点除于总点数就是答案要的概率.
撸起袖子就是干:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | const N = 10000000; for (let i = 0; i < N; ++ i) { let x = Math.random(1.0); 返回[0, 1) 闭包 let y = Math.random(1.0); if ( check(x, y, 0, 0, 0, 1) && check(x, y, 0, 1, 1, 1) && check(x, y, 1, 1, 1, 0) && check(x, y, 1, 0, 0, 0) ) { valid ++; } } console.log(valid); console.log(N); console.log(Math.round(valid / N * 1000) / 1000); |
const N = 10000000; for (let i = 0; i < N; ++ i) { let x = Math.random(1.0); 返回[0, 1) 闭包 let y = Math.random(1.0); if ( check(x, y, 0, 0, 0, 1) && check(x, y, 0, 1, 1, 1) && check(x, y, 1, 1, 1, 0) && check(x, y, 1, 0, 0, 0) ) { valid ++; } } console.log(valid); console.log(N); console.log(Math.round(valid / N * 1000) / 1000);
注意: Math.random(1.0) 返回[0, 1)闭包区间的随机值, 我们不是说过不能等于0嘛, 不能在正方形边上, 不过没关系, 这个误差可以忽略不计, 因为本身这个模型就是随机采样的, 本身就有少量误差!
跑几次的结果虽然都是 0.213 但是前面有效的点数会有小范围浮动, 比如:
2126142
2127333
2125929
....
如果你的N少写几个零, 那么精度就很有可能不够了, 可能会输出 0.212, 0.214 这样很接近的值. 这种方案的时间复杂度是 O(N) , N也就是你的采样点数, 可以根据需求调整, 比如我只想大概知道概率, 不需要很精确, 那么我的N可以是1e4, 这样也大概给出了 0.200, 0.215, 0.216, 0.204 这样的概率 (20%左右), 而实现起来相当的快, 运行速度也是秒出结果. 有它的优势.
当然, 如果你硬是说获取精度若是太高则需要的N很大, 速度很慢, 你是对的, 但我们可以用并行计算来完成, 比如 C++里, 我们可以这么写:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | int main() { srand(time(NULL)); // seed const int N1 = 1000; const int N2 = 100000; int n = 0; int c = 0; Concurrency::critical_section cs; Concurrency::parallel_for(0, N1, [&](int i) { // 并行 int t = check(N2); // 并行小块分别统计 cs.lock(); // 避免同时更新这个变量 n += N2; // 总样本 c += t; // 有效值 cs.unlock(); }); cout << "概率 ~= " << setprecision(9) << (double)c / n << endl; return 0; } |
int main() { srand(time(NULL)); // seed const int N1 = 1000; const int N2 = 100000; int n = 0; int c = 0; Concurrency::critical_section cs; Concurrency::parallel_for(0, N1, [&](int i) { // 并行 int t = check(N2); // 并行小块分别统计 cs.lock(); // 避免同时更新这个变量 n += N2; // 总样本 c += t; // 有效值 cs.unlock(); }); cout << "概率 ~= " << setprecision(9) << (double)c / n << endl; return 0; }
不止C++, 其它语言 如C#里也有 Parallel For, 思路都是一样的, 这里就不具体说了, 可以看我之前写的相关C++博文(英文): C++ Coding Exercise – Parallel For – Monte Carlo PI Calculation
英文: Monte Carlo solution for Mathematics × Programming Competition #7
loading...
上一篇: STEEM API 系列之获取货币转换
下一篇: 渐行渐远的孩子