警告：我认为本文中的内容本身并没有直接的实际应用（除非您是无核战争幸存者，并且需要从头开始构建数学模型或类似的东西）。虽然有时候我喜欢回到基础。这里是π\ pi和没有任何演算或先验功能的弯曲区域的外观。
该算法从简单的数字理论思维开始。一些整数形成纯毕达哥拉斯三元组（x，y，z）（x，y，z），其中x 2 + y 2 = z 2 x ^ 2 + y ^ 2 = z ^ 2。例如3 2 + 4 2 = 5 2 3 ^ 2 + 4 ^ 2 = 5 ^ 2。通过强力搜索很容易找到x 2 + y 2 = 5 2 x ^ 2 + y ^ 2 = 5 ^ 2的所有解，因为我们知道xx和yy不能大于55。在这里，它们是：
0 2 + 5 2 = 5 2 3 2 + 4 2 = 5 2 4 2 + 3 2 = 5 2 5 2 + 0 2 = 5 2 \开始{对齐} 0 ^ 2 + 5 ^ 2＆amp; = 5 ^ 2 \\ 3 ^ 2 + 4 ^ 2＆amp; = 5 ^ 2 \\ 4 ^ 2 + 3 ^ 2＆amp; = 5 ^ 2 \\ 5 ^ 2 + 0 ^ 2＆amp; = 5 ^ 2 \ end {aligned} （加上所有负数组合，但让我们坚持使用非负整数，只计算4个解。）如果我们放宽方程，对x 2 + y 2≤5 2 x ^ 2 + y ^ 2 \在5 ^ 2中，答案为26。为什么要关心？好吧，如果t t是x 2 + y 2≤n 2 x ^ 2 + y ^ 2 \ le n ^ 2的解的总数，则
lim n＆amp; rightarrow∞4 t（n + 1）2 =π\ lim_ {n \ to \ infinity} \ frac {4t} {（n + 1）^ 2} = \ pi或者，在代码中，这很简单估计π\ pi的程序，对于n变量的较大值，其精度更高：
进口标准; ulong sq（ulong x）pure {return x * x; } void main（string [] args）{const n = args。长度＆gt; 1个args [1]。至 ！乌龙：20;乌龙共; foreach（x; 0. n + 1）{foreach（y; 0. n + 1）{如果（sq（x）+ sq（y）＆lt; = sq（n））总计++; }} / * //或者，对于函数式编程爱好者：const total = cartesianProduct（iota（n + 1），iota（n + 1））.count！（p =＆gt; sq（p [0]）+ sq（ p [1]）＆lt; = sq（n））; * / writef（＆＃34;％。12f \ n＆＃34;，4.0 *总计/ sq（n +1））; }
好的，这比22 7 \ frac {22} {7}更准确。不过，与大多数π\ pi公式不同，有一个简单的图表显示了它的工作原理。想象一下，我们以明显的方式在二维网格上布置（x，y）（x，y）整数对（其中x x和y y范围从0 0到n n）。下图显示了n = 10 n = 10的示例，箭头r r从原点指向（6，8）（6，8）。 r r以及x x和y y组成一个直角三角形，因此毕达哥拉斯定理说x 2 + y 2 = r 2 x ^ 2 + y ^ 2 = r ^ 2。对于（6，8）（6，8），r = 10 = nr = 10 = n，因此（6，8）（6，8）在边界上作为x 2 + y 2≤10 2 x的解^ 2 + y ^ 2 \ le 10 ^ 2。该边界（一组实值点与原点的恒定距离n = 10 n = 10）形成了四分之一圆。
圆是简单的凸形，并且网格点的间距均匀，因此四分之一圆内的点数将与面积大致成比例。更具体地说，四分之一圆内所有网格点的分数将大致为四分之一圆的面积除以所有点周围的正方形面积。四分之一圆的面积是πr 2＆amp; div 4 \ pi r ^ 2 \ div 4，在面积r 2 r ^ 2的平方内（请记住，n = rn = r），所以π4 \ frac {\ pi} {所有点中的4}代表解。 x x和y y值从0 0到n n计数，因此有（n +1）2（n + 1）^ 2个网格点。重新排列方程式，您将获得一个用于根据解数估算π\ pi的公式。随着n n的增大，网格点将不断绘制任意精确的圆（就像高分辨率的计算机监视器一样），因此估算值在极限范围内是精确的。
上面的代码很简单但是很慢，因为它蛮力扫描了所有可能的x x和y y值（n + 1）×（n + 1）（n + 1）\ x（n + 1）。但是显然我们不需要扫描所有值。如果我们知道x 2 + y 2≤n 2 x ^ 2 + y ^ 2 \ le n ^ 2，那么减小x x或y y只会给我们另一个解决方案。找到解决方案后，我们无需继续测试较小的值。最终，我们只需要找到边界周围的积分点即可。这是一个基于该想法的更快算法。
想象我们沿着积分x x值扫描，找到最大积分y y值仍然可以为我们提供解决方案。这使我们在下图中用红色标记了边界线。如果对于给定的xx值y = 8 y = 8，我们立即知道给定的xx值有8 + 1 = 9 8 +1 = 9个解（+1 + 1来计算y = 0 y = 0的解） 。
请注意，当x x从0 0扫描到n n时，y y从n n开始并减小到00。重要的是，它仅减小-单调。因此，如果我们从0 0到n n扫描x x，则可以从前一个边界点开始向下搜索来找到下一个边界y y点。这是一些代码：
ulong y = n，总计； foreach（x; 0.。n +1）{while（sq（x）+ sq（y）＆gt; sq（n））y-;总计+ = y + 1; }
此版本仍具有嵌套循环，因此看起来仍然是O（n 2）O（n ^ 2）。但是，内部while循环对每个x x值执行不同的次数。通常y根本不会触发。实际上，由于y从n开始并且单调减小到0，我们知道y--总共将被精确执行n次。该代码中没有指令可以执行超过O（n）O（n）次，因此整个算法为O（n）O（n）。
使用64b ulong整数，在溢出前有效的n的最大值为4294967294：
有一些方法可以使用数值积分技巧来加快收敛速度，但是我喜欢这种算法仅使用整数算术（直到最后除法）的方式，并且可以从简单的图中直接理解。
也许您会觉得有点受骗，因为该算法假设圆面积的πr 2 \ pi r ^ 2公式。当然，这可以说是包含在“高中数学”中的，但是除非学生学习积分并以这种方式推导，否则有些事情会被告知要记住。但是，如果我们要假设πr 2 \ pi r ^ 2，为什么不也假设三角函数理论，而只使用π= 4 atan（1）\ pi = 4 \ atan（1）？
伟大的古希腊数学家阿基米德（Archimedes）在2000年前就算出了没有整数演算（或与此有关的三角函数）的圆形区域。他从对规则（即等边）多边形的精妙见解开始。
欧几里得的《几何元素》（很容易从平行四边形的定理推导出来）的第一本书中就已经有一个著名的证明，即三角形区域的“半基本乘以高度”公式。方便地，任何规则的多边形都可以分成与中心相连的相等三角形。例如，正六边形会分成六个三角形，如下图所示。我们可以采用任意一个三角形（它们都是相同的），并将“底”称为多边形的一面。然后，“高度”是从底部中心到多边形中心的线。
现在，这是阿基米德的巧妙见解：三角形面积与底数的比率为h 2 \ frac {h} {2}。如果将所有面积加起来，就可以得到多边形的面积。同样，如果将所有基准相加，将得到多边形的周长。因为对于所有三角形，三角形的面积/底比都是恒定的h 2 \ frac {h} {2}，所以整个多边形的面积/周长比是相同的h 2 \ frac {h} {2}。根据公式，任何规则多边形的面积为P×h 2 P \ times \ frac {h} {2}（其中P P为周长）。
如果您将圆视为具有无限多个边的规则多边形（以使hh成为圆的半径），并使用圆周长（2πr 2 \ pi r）作为π\ pi的基本定义，则这意味着一个圆的面积是2πr×r 2 =πr 2 2 \ pi r \ times \ frac {r} {2} = \ pi r ^ 2。
当然，阿基米德是一位受人尊敬的数学家，他仅凭假设多边形的真假性就等于一个圆形就无法逃脱（反例：所有多边形的拐角都颠簸，而圆形却不然）。他使用当时流行的那种矛盾的几何证明。 （他甚至进一步研究了球体，圆柱体，抛物线和其他弯曲的物体，几乎在几千年前就发明了像现代实物分析之类的东西。）可悲的是，他的所有数学工作都没有幸存，但是他的测量的关键部分一个圈子。
这是高级版本。阿基米德声称，圆的面积为2πr×r 2 2 \ pi r \ times \ frac {r} {2}。假设您认为他的值太小，圆确实大于2πr×r 2 2 \ pi r \ times \ frac {r} {2}。这意味着圆内有足够的空间容纳正则多边形，该正多边形也大于2πr×r 2 2 \ pi r \ times \ frac {r} {2}。但是阿基米德说这是矛盾的，因为对于任何这样的多边形，h < r h＆lt r，P＆lt; 2πr P \ lt 2 \ pi r（因为多边形的每一边都是一条比连接相同点的圆弧短的直线），所以面积A = P×h 2 < 2πr×r 2 A = P \ times \ frac {h} {2} \ lt 2 \ pi r \ times \ frac {r} {2}。多边形的面积不能大于和小于2πr×r 2 2 \ pi r \ times \ frac {r} {2}。
阿基米德认为，如果您认为2πr×r 2 2 \ pi r \ times \ frac {r} {2}太大，并且圆的面积小于该范围，则存在类似的矛盾。 在这种情况下，他可以制作一个也小于2πr×r 2 2 \ pi r \ times \ frac {r} {2}的多边形，但仍然环绕圆。 对于此多边形，h = rh = r，但他说多边形的周长必须大于2πr 2 \ pi r 1，因此多边形的面积必须大于2πr×r 2 2 \ pi r \ 倍\ frac {r} {2}，即使它也应该更小。 如果这两种情况都导致矛盾，那么我们剩下的唯一选择就是圆圈面积为πr 2 \ pi r ^ 2。