任意の円の内部にある格子点を求めるＣ言語のプログラムコードを以下に示す。
#include <stdio.h>
#include <math.h>
/*天井と床を求める*/
void cf(double x, double r, int *low, int *high){
*low = ceil(x-r);
*high = floor(x+r);
}
int main(void){
double a,b,r,p;
int start,end,top,bottom,i,j,num;
num=0;
printf("a b r =");
scanf("%lf %lf %lf",&a,&b,&r);
cf(a,r,&start,&end);
for(i=start;i<=end;i++){
p = sqrt(pow(r,2)-pow((a-(double)i),2));
cf(b,p,&bottom,&top);
for(j=bottom;j<=top;j++){
num = num + 1;
}
}
printf("中心:(%lf,%lf) 半径:%lf 格子点の数: %d \n",a,b,r,num);
return 0;
}
aは円の中心のx座標,bは円の中心のy座標,rは円の半径を表し、その円の内部にある格子点の数が出力される。使用するときはターミナルからa,b,rを打ち込んでから使う。
以下、このプログラムの説明を書いていく。
まず、最初に定義された関数cfは天井(ceiling)と床(floor)を求める関数である。これだけでは何をいっているわかりにくいので図を用いて説明する。
天井とはある数を超える最小の整数のことである。例えば3.5の天井は4、-2.1の天井は-2である。床とはある数を超えない最大の整数のことである。例えば3.5の床は3、-2.1の床は-3である。
下の図を見てみるとa-rの天井(start)が格子点のx座標の最小値となっていることがわかる。この図の例ではそれは-3である。
また、a+rの床(end)が格子点のx座標の最大値となっていることがわかる。この図の例ではそれは0である。
格子点のx座標は$\rm start\leqq i \leqq end$を満たす整数iに限られることがわかったので、各iごとに格子点の数を求めて足し合わせればよいとわかる。
下の図で三平方の定理より
$$p=\sqrt{r^2-(a-i)^2}$$
となるとわかるので、b+pの床(top)がx座標がiの格子点のy座標の最大値であり、b-pの天井がx座標がiの格子点のy座標の最小値であるとわかる。この図の例ではi=-2のときtop=4,bottom=-1とわかる。
格子点のx座標がiのとき、そのy座標jは$\rm bottom \leqq j \leqq top$を満たすとわかるので各iごとにjの数をカウントアップし、すべてのiについてそれを行えば格子点の数(num)が求められる。
実行結果は以下の通り
中心:(3.000000,2.000000) 半径:6.000000 格子点の数: 113
中心:(3.450000,-9.310000) 半径:5.510000 格子点の数: 93
中心:(0.000000,0.000000) 半径:1.100000 格子点の数: 5
備考
半径rの円内部の格子点の数をN(r)とすると
$$\lim_{r \to \infty}\frac{N(r)}{r^2}=\pi$$
という面白い関係式が成り立つそうなので$r=10^4$としてみる。実行結果は
中心:(0.000000,0.000000) 半径:1.000000e+004 格子点の数/半径の2乗:3.141591
となり、確かに$\pi$に収束しそうだとわかる。