임의의 원의 내부의 격자점의 수를 구한다

임의의 원의 내부에 있는 격자점을 구하는 C 언어의 프로그램 코드를 이하에 나타낸다.

#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$에 수렴할 것 같다는 것을 알 수 있다.