prog/02 のバックアップ(No.3)

バックアップ一覧
差分を表示
現在との差分を表示
ソースを表示
prog/02 へ行く。
- 1 (2015-11-11 (水) 14:58:21)
- 2 (2015-11-11 (水) 17:31:11)
- 3 (2015-11-13 (金) 20:32:27)
- 4 (2015-11-13 (金) 20:32:27)
- 5 (2015-11-21 (土) 13:42:20)

円周率の計算プログラムはどこまで小さくなるか

(by K, 2015.11.11)

↑

きっかけ

これはprog/01の続きです。
hikaliumさんがprog/01を見て、なんとtwitterに投稿できるレベルにまで文字数を減らしてしまいました。
- https://twitter.com/hikalium/status/664052511121145856
```
#include<stdio.h>
int main(){
int N=20000,i=N-1,j,c=0;
for(j=1;j<N;j++){
while(i*i+j*j>N*N)i--;c+=i;
}
printf("pi=%d\n",c+N-1);
return 0;
}
```
おおこれはすごい！そしてこれを見たら私ももっとやってみたくなってしまったのです（笑）。
- もはや当初の目的の「円周率の計算方法を中学生に分かりやすく説明する」を激しく逸脱しています。

↑

スタートはここから

上記のhikariumさんのコードをスタートにしてもいいのですが、それだとそこまでの過程がすぐには分からないかもしれないので、ここはあえてprog/01の終わりから始めることにします。
```
#include <stdio.h>

#define N 10000

int main()
{
    int i = N - 1, j, c = 0;
    for (j = 1; j < N; j++) {
        while (i * i + j * j > N * N) i--;
        c += i;
    }
    printf("pi=%d\n", c * 4 + 4 * N - 2);
    return 0;
}
```

↑

4倍の計算を消す

prog/01の「整数演算だけにする」に書いた通り、cはN*N倍で計算されます。ということは、Nを2万にしてやれば、cの値はさらに4倍で計算されて、結果的にprintfでの4倍処理を省略できます（hikaliumさんのアイデアです）。・・・そしてこのとき、-2の計算はむしろやらないほうがパイの値に近づくので、これをやらないことにします。

またその他の些末な変更も加えて、とにかくプログラムを小さくしてみます。

#include <stdio.h>
main()
{
    int N = 20000, i = N - 1, j, c = 0;
    for (j = 1; j < N; j++) {
        while (i * i + j * j > N * N) i--;
        c += i;
    }
    printf("%d\n", c + N);
}

↑

初期値を工夫する

最後にcに対してNを加算していますが、だったら最初のc=0をc=Nにすればいいと思いました。

あとiの初期値はN-1になっていますが、仮にここがNだったとしても、すぐにwhileループでN-1に修正されるので問題ありません。

#include <stdio.h>
main()
{
    int N = 20000, i = N, j, c = N;
    for (j = 1; j < N; j++) {
        while (i * i + j * j > N * N) i--;
        c += i;
    }
    printf("%d\n", c);
}

↑

ライバル登場

ここで、ライバルを登場させます（笑）。短いコードで円周率を計算すると言ったら、定番はこれだと思うからです。
```
#include <stdio.h>
main()
{
    int i;
    double p = 0;
    for (i = 1; i < 4000000; i += 4)
        p += 4. / i - 4. / (i + 2);
    printf("%.8f\n", p);
}
```
うーん、やっぱり短いですね。これは三角関数のタンジェントの逆関数の多項式展開の公式を使っているのですが、こんなの中学生向けには説明できません。だからmath/03では紹介しませんでした。
上記の400万は、N=2万の結果とほぼ同じくらいの精度を出すために必要なもので、つまりこれは結構重い計算になります。大雑把にいって、10倍くらい遅いです。
上記の計算では倍精度の浮動小数点を使っていますが、それは必須ではありません。pを1億倍とかで計算することにして式を整理すれば、整数演算だけでも同等の計算はできます。
しかし割り算はなくせません。割り算は四則演算の中では一番難しいものなので、それすら使わないprog/01のアルゴリズムは優秀だと思います。

うーん、それにしても文字数で負けているのは悔しい。mainの中身がどのくらい違うのか比べてみよう。

int N=20000,i=N,j,c=N;for(j=1;j<N;j++){while(i*i+j*j>N*N)i--;c+=i;}printf("%d\n",c);
int i;double p=0;for(i=1;i<4000000;i+=4)p+=4./i-4./(i+2);printf("%.8f\n",p);

おお、8文字しかちがっていない！でも76文字に対しての8文字だから、それでも1割以上負けていることになるけれど・・・。

四則演算以上のものを使ってもよければ、<math.h>をincludeした上で printf("%.8f\n", atan(1)*4); という面白くもなんともない方法で円周率は計算できます。というかそもそもx86のCPUは定数パイを求める命令があるので、それを使えば1命令です（単に記憶しているのを読み込んでいるだけだけど）。・・・でも、僕たちがやりたいことってそういうことじゃないですよね。だからそういう無粋な方法は検討すらしないことにします。

↑

そしてOSECPU-VMへ

ここまできたら、OSECPU-VMで限界を見てみたくなってしまいました。割り算を使わないバージョンは、最終的に39バイトまで行きました。

#include "osecpu_ask.h"

#define N   20000

    Int32s n:R00, i:R01, c:R02, nn:R03, j:R04;
    n = N;
    i = n;
    c = i;
    nn = c * c;
    for (j = 1; j != n; j++) {
        for (;;) {
            if (j * j + i * i <= nn) break;
            i--;
        }
        c += i;
    }
    api_putStringDec('\1', c, 9, 0);

一方の割り算利用バージョンは、以下のようになってなんと31バイトでした。
Nがやたらと大きくなったのは、整数演算版だとこれくらいやらないと上記と同じくらいの精度が出ないということと、2のベキ定数だとアプリのバイナリサイズが小さくなるからです。
- 最初は浮動小数点バージョンも試作しましたが、努力しても40バイトすら切れず、これは見込みがないとさっさと整数版に取り組みました。
```
#include "osecpu_ask.h"

#define N   0x8000000

    Int32s c:R00, i:R01, p:R02;

    c = 400000000;
    for (i = 1; i < N; i += 2) {
        p += c / i;
        c *= -1;
    }
    api_putStringDec('\1', p, 9, 0);
```
(2015.11.13追記) 31バイトのバージョンをベースにさらにVM側の仕様も改良したところ、ついに21バイトになってしまいました。

↑

こめんと欄

コメント	お名前	NameLink