这道小学六年级的数学题，恕我直言没几个人会做

godweiyang

发布于 2021-12-21 20:33:48

4860

发布于 2021-12-21 20:33:48

文章被收录于专栏：算法码上来

作者 | godweiyang

出品 | 公众号：算法码上来（ID：GodNLP）

今天网上冲浪的时候突然看到一道小学六年级的数学题，如上图所示，求阴影部分的面积。

我下意识就想到了微积分，这不就建立坐标系，求出交点，计算积分就行了嘛。转念一想，小学生哪里会积分，这道题一定有简单的解法。

经过一番努力，我尝试出了一共两种精确解法和一种近似解法，但都超出了小学生的理解范畴。

首先我们把图左右翻转一下，方便建立坐标系，如下图所示：

积分方法

首先联立半圆和对角线的方程：

\begin{aligned} x^2+(y-4)^2-16&=0 \\\ x+2y-4&=0 \end{aligned}

解出交点E的坐标是

(\frac{12}{5}, \frac{4}{5})

所以可以直接求下面积分，得到要求的阴影面积：

\int_{0}^{\frac{12}{5}}{\left(\sqrt{16-x^2}+4\right)}dx + \int_{\frac{12}{5}}^{4}{\left(2-\frac{1}{2}x\right)}dx

利用简单的高等数学知识，可以解出答案是下面这样的，约等于1.252。

\frac{32}{5} - 8 \arcsin{\frac{3}{5}}

但是小学生肯定不会积分呀，就算高中生大多数都不会！

间接计算

既然我们不方便算这个不规则图形的面积，我们可以计算其他规则图形的面积，然后加减之后间接得到答案。

图中阴影部分面积其实等价于三角形ACE+梯形BOCE-扇形BOE，而这三部分其实都非常好算。

因为E的坐标是

(\frac{12}{5}, \frac{4}{5})

，所以最难算的就是扇形BOE的面积，也就是要求角OBE的大小。很轻松就可以得到角度是

\arcsin{\frac{3}{5}}

，所以三部分运算后可以得到和上面积分相同的答案。

这个方法初中生都会做了！但是对于小学生来说难度还是有点大了。

蒙特卡洛采样方法

如果往第一象限正方形内随机抛撒若干豆子，然后看掉落在阴影部分里的豆子所占的比例，那么就可以近似估计出阴影部分的面积了，这就叫蒙特卡洛采样。

那我们可以随机生成两个随机数，作为采样点的坐标，然后判断是否在阴影部分里，重复上述过程，记录次数就行了。

我这里写了C++和CUDA两个实现，CUDA可以并行采样，速度快很多。

执行的话，C++保存为test.cpp，然后执行g++ test.cpp -o test和./test。CUDA保存为test.cu，然后执行nvcc test.cu -o test和./test。

实测CUDA相同时间内可以采样的次数大约是100倍左右，主要受到原子写操作的局限，不然可以达到1000倍左右的加速。

但是小学生会CUDA吗？显然不可能。

C++实现

#include <iostream>

typedef unsigned long long int ull;

int main() {
  srand(time(NULL));
  ull res = 0;
  ull cnt = 1e8;
  for (ull i = 0; i < cnt; ++i) {
    double x = double(rand()) / RAND_MAX;
    double y = double(rand()) / RAND_MAX;
    if (x * x + (y - 1) * (y - 1) > 1 && x + 2 * y - 1 < 0) {
      res += 1;
    }
  }
  printf("%.5f\n", res * 16.0 / cnt);
  return 0;
}

CUDA实现

#include <iostream>
#include <chrono>
#include <curand_kernel.h>

typedef unsigned long long int ull;

__global__ void monte_carlo_kernel(ull* res, ull cnt, const int seed) {
  ull i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i >= cnt) return;
  curandStatePhilox4_32_10_t state;
  curand_init(seed, i, 0, &state);
  double2 rand = curand_uniform2_double(&state);
  double x = rand.x;
  double y = rand.y;
  if (x * x + (y - 1) * (y - 1) > 1 && x + 2 * y - 1 < 0) {
    atomicAdd(res, 1);
  }
}

void launch_monte_carlo(ull* res, ull cnt, cudaStream_t& stream) {
  ull grid_dim = cnt >> 10;
  monte_carlo_kernel<<<grid_dim + 1, 1024, 0, stream>>>(
      res, cnt,
      std::chrono::duration_cast<std::chrono::microseconds>(
          std::chrono::system_clock::now().time_since_epoch())
          .count());
}

int main() {
  ull* res;
  cudaMallocManaged(&res, sizeof(ull));
  res[0] = 0;
  ull cnt = 1e10;
  cudaStream_t stream;
  cudaStreamCreate(&stream);
  launch_monte_carlo(res, cnt, stream);
  cudaDeviceSynchronize();
  printf("%.5f\n", res[0] * 16.0 / cnt);
  cudaFree(res);
  return 0;
}