wayne 发表于 2017-9-24 23:59
先占个坑。
对于给定分辨率的“孔方兄”，一个像素格子可以由其中心点的坐标以及边长决定。标记成三元组{a ...
计算出来每迭代一次，小正方形的个数乘以8,最终格子中心落在圆内的个数是{{1,0},{2,12},{3,136},{4,1056},{5,8536},{6,68488},{7,547928}}
除以对应的总数，得到圆内的格子的占比比例，{{1,0},{2,3/16},{3,17/64},{4,33/128},{5,1067/4096},{6,8561/32768},{7,68491/262144}}
数值结果就是{{1.,0.},{2.,0.1875},{3.,0.265625},{4.,0.257813},{5.,0.260498},{6.,0.261261},{7.,0.261272}}
极限值看上去像是$ 0.261272$
顺便贴一下我画的图
{3,17/64},{4,33/128},
显然不对，圆内的应该只增不减
另外最好再给出圆上的数目，这样就可以非常容易的同时算出上下界
mathe 发表于 2017-9-25 16:53
……………………
最好再给出圆上的数目，这样就可以非常容易的同时算出上下界
取上下界的中值，即圆周上的非白方块只计一半，收敛速度会快得多吧，很早就能看到大致结果。
老大们说的对，重新计算了下：
kernel=DeleteCases,{0,0}];
kernel={{-1,-1},{0,-1},{1,-1},{1,0},{1,1},{0,1},{-1,1},{-1,0}};
rect={{-1,-1},{-1,1},{1,1},{1,-1}};
g:=Module[{ps=points},{Flatten]/6 #+k]&/@rect];(Sign]]+Sign]])/2}],{k,ps[]/3 #+p[]&/@kernel}],{p,ps[]}],1],ps[]/3}]
pp=Nest;
Graphics[{EdgeForm],Flatten]==-1,Gray,j[]==0,Blue,j[]==1,Red],Rectangle]+{-1,-1}*pp[]/2,j[]+{1,1}*pp[]/2]},{j,pp[]}],1],Thickness[.005],RGBColor,Circle[{0,0},1/2]}]
统计所有的圆外，圆上，圆内的格子数目分别如下：
{{1},{1.},8}
{{1/16,7/16,1/2},{0.0625,0.4375,0.5},64}
{{13/64,19/128,83/128},{0.203125,0.148438,0.648438},512}
{{123/512,51/1024,727/1024},{0.240234,0.0498047,0.709961},4096}
{{259/1024,145/8192,5975/8192},{0.25293,0.0177002,0.72937},32768}
{{2117/8192,389/65536,48211/65536},{0.258423,0.00593567,0.735641},262144}
{{34119/131072,1045/524288,386767/524288},{0.26030731,0.0019931793,0.73769951},2097152}
{{547243/2097152,2801/4194304,3097017/4194304},{0.26094580,0.00066781044,0.73838639},16777216}
以不同的颜色标记，画图如下：
@wayne 你的结果还是有问题吧。
所谓圆内/外，指的是整个小方块完全在圆内/外，剩下的就是被圆周所穿过的方块。
按｛迭代次数，{圆上，圆内}，{圆内+圆上/2,"面积"}}来陈列，前 3 项数据人工也容易数得出，为
{1, {8, 0}, {4, 1/2}}
{2, {28, 32}, {46, 23/32=0.71875}}
{3, {76, 332}, {370, 185/256=0.72265625}}
第4次迭代
{4, {204,2908}, {3010, 1505/2048=0.73486328125}}
@wayne 你是按小方块的中心来划分的么？
看wayne给出的图，蓝色圆(边界上的)太多了
这道题目还可以用Bresenham算法求出和圆相交各小正方形。而对于每个这样小正方形，实际上可以使用它两坐标的三进制形式计算出其下圆内黑格数目。这样做的好处是空间复杂度比较小
这个题目还可以用蒙特卡洛法，然后我计算划分到20层后的蒙特卡罗法，结果在0.736左右
本帖最后由 lsr314 于 2017-9-26 18:29 编辑
不考虑与圆相交而被切割的误差，只以小正方形的中心位置来划分圆内圆外，那么可以用三进制来判断，代码如下：
f :=
MemberQ[(IntegerDigits, 3] - 1)^2 +
Take, 3] -
1, -Length, 3]]]^2, 0]
ff := (s = t = 0;
Do, ,
If[((i + 1/2)/3^l - 1/2)^2 + ((j + 1/2)/3^l - 1/2)^2 < 1/4,
s = s + 1, t = t + 1]], {j, 0, 3^l - 1}, {i, 0, 3^l - 1}];
s/(s + t))
Do]}], {l, 7}]
结果如下：
{1,1.}
{2,0.8125}
{3,0.734375}
{4,0.742188}
{5,0.739502}
{6,0.738739}
{7,0.738728}
本帖最后由 chyanog 于 2017-9-26 20:35 编辑
写了两个版本的Mathematica代码：
第一个需要11.1版本，使用了MengerMesh和计算几何相关的函数，代码简单易懂，不过效率不高
Clear["`*"];
n=3;
poly=MeshPrimitives,2];
Graphics[{
Antialiasing->False,
EdgeForm,
Red,Select,#]&],
Blue,Select,#]&],
Gray,Select,#]&&RegionDisjoint,#]&],
Thick,Magenta,Circle[{0.5,0.5},0.5]
}]
第二个版本速度快一些，不限于版本11
Clear["`*"];
carpet:=Nest&,1,n];
n=3;
ArrayPlot,Mesh->All]
p=Rescale[{{#1,#2},{#1+1,#2},{#1+1,#2+1},{#1,#2+1}}&@@@N@Position,1]];
Graphics[{
EdgeForm,
Antialiasing->False,
Red,Polygon@Select<0.5&]&],
Blue,Polygon@Select&/@#},Min@m<0.5&&Max@m>0.5]&],
Gray,Polygon@Select>0.5&]&],
Thick,Magenta,Circle[{0.5,0.5},0.5]
}
]
补充内容 (2017-9-27 10:46):
次数，圆上，圆内，圆外，“面积”
{1, 8, 0, 0, 0.5}
{2, 28, 32, 4, 0.71875}
{3, 76, 332, 104, 0.72265625}
{4, 204, 2908, 984, 0.73486328125}
{5, 580, 23900, 8288, 0.73822021484375}