C言語 円内の整数格子点とπ

Question

原点を中心とする半径nの円周または円内にある整数格子点の数を数えます。
そののち、整数格子点の数/(n^2)を計算しπに近づいているかを確かめる、というものです。
一応計算はできますが、計算時間が遅いです。
1億までを3秒で計算しなければいけません。
なにか良い方法があれば、教えてください。

#include<stdio.h>
#include<math.h>

int i=0,j;
long check;
double answer;

double count(int m){
for(i=1;i<m;i++){
for(j=1;j<m;j++){
if(ii+jj<=mm){
check++;
}
}
}
answer=((double)(check4+m4+1)/(double)(mm)); //long型とint型の計算

return answer;

}

int main(void){

int n;
scanf("%d",&n);
printf("%d:%.15f\n",n,count(n));

return 0;

}

Comments

[ikuwow]

解決済みですが、今後見る人のためにもソースコードをMarkdown記法でシンタックスハイライトして見やすくしていただけると助かります。

Answer 1

求めたい格子点のうち、
x座標x=iであるものの個数は、√(m^2 - i^2)の整数部分から簡単に求まります。
あとはi方向でループするだけです。

もしくは図のように走査します。(歪んでますが)

他にもいろいろあると思うので
頑張って頭をひねってください。
方針は
・ループ回数を減らす。
・使いまわせるものは使い回す。
です。

Comments

[marikot]

プログラムはまだ完成していませんが、yahoo知恵袋よりもわかりやすく説明していただいたので、ベストアンサーにさせていただきました。アルゴリズムは理解できたので
実際にコードに起こすのは自分で頑張ってみようとおもいます。

Answer 2

すでに解決済みではありますが、、
とても興味深いテーマでしたので、私なりに考えてみました。

---

まず、以下の図をご覧ください。

青い扇形は半径R（Rは整数）の円の第一象限です。
この図に収まる、辺の長さが1の正方形（以下、単位正方形と記載します）の個数を考えます。

赤い部分（①）は、上の図に収まる、各辺の長さがもっとも大きな整数を持つ正方形です。

②、③、④は、①の正方形の辺を伸ばして第一象限を4分割した、残りの部分です。

このうち、②と③は合同になります。

④に収まる単位正方形は存在しません。

つまり、第一象限に収まる単位正方形の個数は

① + 2 * ②

で求めることができます。

これを4倍すると、半径Rの円に収まる単位正方形の個数となります。

さらに、円に収まる格子点の数は、単位正方形の個数にX軸、Y軸上の格子点の数を足せばよいので、

半径`R`の円に収まる格子点の数 = 4 * (① + 2 * ②) + 4 * R + 1 // 最後の"+ 1"は原点

となります。

①の面積は

floor(R / √2) ^ 2

ですので、これが①に収まる単位正方形の数です。

②についてはozwk様に紹介いただいた方法で求めることができますが、xの値は

floor(R / √2) + 1 <= x < R

の範囲で考えれば十分なので、

② = Σ floor(√(R^2 * x^2)) : floor(R / √2) + 1 <= x < R

と、なります。

まとめると、

半径`R`の円に収まる格子点の数 = 4 * (floor(R / √2) ^ 2 + 2 * (Σ floor(√(R^2 * x^2)))) + 4 * R + 1

です。

式で書くと意味不明ですねw

以下が、上の式をJavaで実装したものです。
（すみません、C言語は不勉強です。。）

public class Main {
    
    private static final long R = 100000000;
    
    private static final long SQUARED_R = R * R;
    private static final long R_DEVIDED_BY_SQRT2 = (long) (R / Math.sqrt(2));
    private static final long LARGEST_SQUARE = R_DEVIDED_BY_SQRT2 * R_DEVIDED_BY_SQRT2; // ①

    /**
     * @param args
     */
    public static void main(String[] args) {
        
        System.out.println("R\t= " + R);
        System.out.println("R * R\t= " + SQUARED_R);
        System.out.println("R / √2\t= " + R_DEVIDED_BY_SQRT2);
        System.out.println("第一象限に収まる最も大きな正方形の面積\t= " + LARGEST_SQUARE);
        System.out.println();
        
        long temp = 0; // ②
        
        for (long x = R_DEVIDED_BY_SQRT2 + 1; x < R; x++) {
            long y = (long) Math.sqrt(SQUARED_R - x * x);
            
            temp += y;
        }
        
        long result = 4 * (2 * temp + LARGEST_SQUARE) + 4 * R + 1;
        
        System.out.println(result);
    }

}

実行結果

R	= 100000000
R * R	= 10000000000000000
R / √2	= 70710678
第一象限に収まる最も大きな正方形の面積	= 4999999983219684

31415926535867961

以下の環境で実行したところ、1秒弱で結果を得ることができました。

Windows7 Home Premium
Intel Core i7-2600 3.40GHz
JDK 1.6

もっとも、丸め誤差などを一切、考慮していないので、
この計算が合っているかはちょっと、自信がありませんが。。

Comments

[marikot]

すごくわかりやすい説明ありがとうございます。
確かにloopの計算が少なくてすむので、確実に早く計算できますね！
floorを使ったことがなかったので、コードを書いて確認してみようと思います！

Answer 3

ozwkさんのように三平方の定理などでループ回数を減らせば現実的な時間でできるかもしれませんが、それをやるなら円の面積の公式を使えばループしないで計算できるということになってしまいます……。

この手の仕事にはモンテカルロ法が真っ先に思い浮かぶのですが、全ての点を1点1点数えることが重要なのであれば、1億×1億のループを現実的な時間で処理するのは不可能ですね。

仮にCPUクロックが10GHzの超高速マシンをお使いになっているとして、ループ1回分に10クロックかかるとして、1億×1億の処理時間を計算すると、

100,000,000 * 100,000,000 / (10,000,000,000 / 10) = 10,000,000 [秒]

約115日と17時間ですね。

追記

ちょっと勘違いしていたみたいです。面目ない。

Comments

[Chironian]

> 円の面積の公式を使えばループしないで計算できる
あっ、元プログラムを見ると、最後に有効桁数15桁に丸めてますね。
ならば、doubleの有効桁数も15桁だし。これが解になるのでは？

それと蛇足ですが、1億は100,000,000ですよ。現実的ではないと言う意味では結果は変わりませんけど。

[catsforepaw]

なんと！　お恥ずかしい。
こっそり書き直しておきます。

[ozwk]

> 円の面積の公式を使えばループしないで計算できる
円の面積の公式で、格子点の個数を算出できますか？

[catsforepaw]

あ、「格子点の数を数える」という分野(?)の話でしたね。勘違いしたかもしれません。結果からπが求まればいいみたいな考えでした。すみません。学が低いもので。
ならば、半径を対角線とする正方形部分を除けば、だいぶループが端折れますね。

[Chironian]

こんばんは。

単位面積に１格子点と考えればできそうな気がしてましたが、矩形の集まりの面積を求める必要があるので誤差が出ます。でも、丸めてるから誤差でてもいいかな？とか思ったのですが、53ビットなら丸めずに済みますね。
残念。

Answer 4

こんにちわ。

直線 y = - x + m と 円 x^2 + y^2 = m^2 で囲まれた領域のみを上記アルゴリズムで
チェックしてみてはいかがでしょうか?

y軸、x軸、直線 y = - x + m で囲まれた領域、すなわち二等辺直角三角形の領域の格子点の数は
上記アルゴリズムを用いるまでもなく簡単に求められますので。

Answer 5

count(), count8(), count_x() の 3 つのメソッドをかいてみました。 (最終敵はコードは count_x() です)
counx_x() は m = 一億で 0.9 秒でした。
実行環境は以下です。
MacBook Pro (13-inch, Mid 2012), OSX 10.11.1, gcc 4.2.1)
(もっと性能のよいマシンなら,さらに短くなるはず。最新の上位機種の MacBook が欲しいなぁ ...)

// See https://teratail.com/questions/20348

#include<stdio.h>
#include<math.h>

double count_x(long m) {
  long count = 0;
  long mm = m * m;
  long i = 0;
  long j = 0;
  long out_j = m;
  for (i = 1; i < m; i++) {
    for (j = out_j; i < j; j--) {
      if (i * i + j * j <= mm) {
        count = count + j - i;
        out_j = j + 1;
        break;
      }
    }
  }

  long d = 0;
  for (d = m - 1; d > 0; d--) {
    if (d * d * 2 < mm) {
      break;
    }
  }
  return ((double)(count * 8 + m * 4 + d * 4 + 1) / (double)(mm));
}

double count(long m) {
  long count = 0;
  long mm = m * m;
  for (long i = 1; i < m; i++) {
    for (long j = 1; j < m; j++) {
      if (i * i + j * j <= mm) {
        count++;
      }
    }
  }
  return ((double)(count * 4 + m * 4  + 1) / (double)(mm));
}

double count8(long m) {
  long count = 0;
  long mm = m * m;
  long i = 0;
  for (i = 1; i < m; i++) {
    for (long j = i + 1; j < m; j++) {
      if (i * i + j * j <= mm) {
        count++;
      }
    }
  }

  long d = 0;
  for (d = m - 1; d > 0; d--) {
    if (d * d * 2 < mm) {
      break;
    }
  }
  return ((double)(count * 8 + m * 4 + d * 4 + 1) / (double)(mm));
}

int main(void) {
  // int n;
  //scanf("%d",&n);
  long n = 10000 * 10000 * 4;

  //printf("%ld: %.15f\n", n, count(n));
  //printf("%ld: %.15f\n", n, count8(n));
  printf("%ld: %.15f\n", n, count_x(n));
  return 0;
}

それぞれの書き換え経緯を書いていきます。

1. count()

==============
まず、質問文のコードを int でなく、long で書き換えました。
このコードは 円の 1/4 のついてだけ計算しています。
（ $ rime ./a.out として実行)

m = 100000, 200000 での結果
100000: 3.141592545700000
real 0m28.971s
200000: 3.141592610525000
real 1m55.339s

m が 2 倍になったとき、処理時間が 4 倍となっています。
このことは後で対処することにします。

1. count8()

=============
まずは、円の 1/8 についてだけ計算するように書き換えました。
（これで処理時間は半分になるはずです)

m = 100000, 200000 での結果
00000: 3.141592545700000
real 0m14.891s
200000: 3.141592610525000
real 1m0.558s

処理時間は半分になりました。

1. count_x()

=============
次は m が 2 倍になったとき、処理時間が 4 倍になっていることに対処します。
処理時間が 2 倍に抑えられていないのは、j のループが無駄なのです。
i のループで 円周上を左から右に移動していきます。2 つのとなりあった 縦線上の円周上の j は
同じか、-1 変化するだけのはずです。
したがって、j を見つけのは、2 回だけの計算で済むのです。
(左隣りの j の値か、 j -1 だけが候補であり、大量にループを回すのは無駄)
この方針で作成したものの処理時間を計測します。

100000: 3.141592545700000
real 0m0.008s

時間が短くなって計測誤差程度になってしまったので、m を増やします。
10000000: 3.141592653505890
real 0m0.128s
20000000: 3.141592653565913
real 0m0.209s
m が 2 倍になったとき、処理時間も 2 倍程度におさえられました。

いよいよ m = 1 億で計測します。
100000000: 3.141592653586796
real 0m0.861s
200000000: 3.141592653589560
real 0m1.709s
400000000: 3.141592653589088
real 0m3.406s

Answer 6

パイの近似値を求めるアルゴリズムなんですね。
Yahoo知恵袋？に同じような質問がありました。
解き方もここのベストアンサーの方法でいいのではないでしょうか？

Answer 7

何点か確認したいことと、アドバイスがあります。

1. 「1億まで」というのは、1辺が1億で間違いないでしょうか（1万×1万で点のほうが1億、ということではないですよね？）。
2. 1億×1億の1兆個の点をまじめに調べていては、3秒では終わりません。うまく工夫して、全点を1つ1つ調べずに数える必要があります。
3. コンピューターの環境にもよりますが、1億×1億=1兆は、int型に入りきらないことがあります。環境によってはlongに入りますが、そうでない場合はさらに対策が必要です。

Comments

[marikot]

説明不足申し訳ありません。
1辺が1億までです。