小赖子的英国生活和资讯

数学 × 程式编写比赛 (第七回) 蒙特卡罗撒点

JustYY.com 小赖子的英国生活和资讯
阅读 桌面完整版
math-and-programming 数学 × 程式编写比赛 (第七回) 蒙特卡罗撒点 I.T. SteemIt 数学 程序设计

math-and-programming

在一个正方形内随机选取一点, 并将此点与四个顶点连上直线, 从而将正方形分割为四个三角形. 求四个三角形之中所有内角均不超过 120° 的概率, 答案准确至小数点后 3 位.

@tvb 在她的方案中 讲了的数学方法需要用到各种三角公式来计算面积, 我数学功底不行, 这种方法对我来说太麻烦(虽然算面积这种想法我也想到了).

暴力, 我想大家都会首先想到, 而且程序好写, 需要注意的一点就是你起始点不能在正方形边上, 然后 x 和 y 坐标 都以一个非常小的 步长 增量, 比如 x += 1e-5, y += 1e-5.

补一句: 这里的暴力并不是真的暴力, 因为你无法穷举所有的点, 这里只是尽可能的穷举“所有”, 有一点小误差.

这样一来, 程序的思路就非常明确了, 只需要判断一下是否是有效点, 然后统计一下除于总个数就是概率.

判断逻辑也较为简单, 根据 (x, y) 连接四个点, 然后得到4个三角形, 依次对每个三角形进行判断内角, 这里我们需要用到一个三角公式.

triangle

已知 三边长, 可以通过 公式: 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

强烈推荐

微信公众号: 小赖子的英国生活和资讯 JustYYUK

阅读 桌面完整版

相关文章：

  1. 你给SteemIt中文微信群拖后腿了么? 阅读 桌面完整版 这年头不缺算法, 就缺数据. 这两天花了很多时间在整API上, 整完之后自己用了一下还觉得真是挺方便的. 今天就突然想看一看自己是否给大家拖后腿了, 于是调用每日中文区微信群排行榜单的API, 刷刷拿着 NodeJs 练手: 1 2 3 4...
  2. Microbit 编程: 简易人工智能让电脑玩游戏 阅读 桌面完整版 上一周, 我们在 Microbit 上编写了第一个游戏: 吃苹果. 我两儿子很喜欢玩, 他们互相比着最高分, 大概极限是35分. 吃苹果的游戏代码和Microbit模拟器: https://makecode.microbit.org/_DV93uT7i0WuK 可能有极限吗? 即使我们做出足够迅速的反应并且没有犯错, 会出现怎么样也无法抓住苹果的情况吗?...
  3. Javascript 中 sleep 函数实现 阅读 桌面完整版 Javascript 中并没有 built-in 的 sleep 函数支持, 在 async/await/Promise 的支持之前, 我们可以用 busy-waiting 的方式来模拟: 1...
  4. 如何在海外 在线观看 2016 欧洲杯比赛? 阅读 桌面完整版 四年一度的欧洲杯足球又开始了, 夏天来了 没有什么比边看足球边喝啤酒更爽的事了. 通过以下链接可以 免费在线观看 2016 欧洲杯比赛 (英语解说). 点击以下链接能在线免费看 2016 欧洲杯比赛 (The following...
  5. 《Steem 指南》之 justyy 在线工具与 API 系列 – 同时给多个帐号发送SBD或者STEEM 阅读 桌面完整版 同时给多个帐号发送SBD或者STEEM STEEMIT 和 BUSY 的前端都有一个内置的钱包工具, 您可以一次给一个帐号发送 SBD 或者 STEEM. 当我们要给很多很多人发送钱的时候, 就显得有些不方便了. 这时候可以用这个在线工具: https://steemyy.com/wallet-tool/...
  6. 来来来 加我公众号玩 石头剪刀布 阅读 桌面完整版 写着玩的, 不过虽然猜拳的策略是完全随机 但是真的很难上 10分. 赢了得1分, 输了扣1分. 平局分数不变. 输入分数或者score查看当前得分. 这个分数会直接记入终身档案中.. 一个微信帐号一个分数. 暂时不支持 表情符号. 微信服务器...
  7. 今天用了机器人给经理请了假 The Zenefit Bot on Slack 阅读 桌面完整版 I used the Zenefit Bot on Slack today to apply for a day...
  8. 为了长远考虑, CN 低保银行降息计划 阅读 桌面完整版 办银行一个月共发了75 SBD 的利息后, 我在想, 这事完全就是赔钱买吆喝. 我们按 7000 SP 来算, 每天利息费用是2.8 SBD, 这么多SP产生的价值(按100%)来算, 一次是(理想状态)...
Exit mobile version