古典的レイトレーサーの実装(第3回)
複数の物体の描画, 影の計算, SSAA, 並列化, AOを実装する

今日やること

  • 複数の物体を扱えるようにする
  • 影の計算
  • 乱数の生成
  • スーパーサンプリング(SSAA)
  • 並列化計算
  • アンビエントオクルージョン(AO)

複数の物体を扱えるようにする

前回までで、一つの球に影をつけて表示することができるようになりました。今回は複数の球を扱えるように拡張してみましょう。

Accelクラス

球の集合を保持するAccelクラスを作りましょう。

メンバ変数

  • std::vector<std::shared_ptr<Sphere>> spheres 球の集合

メンバ関数

  • void add(std::shared_ptr<Sphere> p) 球を追加する
  • bool intersect(const Ray& ray, Hit& hit) const rayと衝突計算を行い、結果をhitに格納する
#ifndef ACCEL_H
#define ACCEL_H
#include <vector>
#include <memory>
#include "ray.h"
#include "hit.h"
#include "sphere.h"
class Accel {
    public:
        std::vector<std::shared_ptr<Sphere>> spheres;

        Accel() {};

        void add(std::shared_ptr<Sphere> p) {
          ...
        };

        bool intersect(const Ray& ray, Hit& hit) const {
          ...
        };
};
#endif

intersectの実装

まず全物体に対して衝突計算します。

Hit hit_each;
for(auto p : spheres) {
  if(p->intersect(ray, hit_each)) {
    //当たった
  }
}

各物体の衝突情報はhit_eachに格納されます。

複数の物体との衝突計算を行う場合、返す必要のある衝突情報は 最も衝突距離が近い物体の衝突情報 だけです。

衝突情報の計算
衝突情報の計算

これを踏まえるとAccel::intersect()は次のようになるでしょう。

bool intersect(const Ray& ray, Hit& hit) const {
    bool isHit = false;

    //hitの距離を大きな値で初期化しておく
    hit.t = 10000000;

    //全ての物体に対し衝突計算をする
    Hit hit_each;
    for(auto p : spheres) {
        if(p->intersect(ray, hit_each)) {
            isHit = true;
            //衝突距離が短いものが見つかったらhitをそれで置き換える
            if(hit_each.t < hit.t) {
                hit = hit_each;
            }
        }
    }

    return isHit;
};

最初にhit.tに大きな値をセットしておき、もしhit_each.thit.tより小さかったらhithit_eachで更新するという処理です。

以上のことを踏まえてAccelクラスを自分で実装してみましょう。

練習問題

  • Accelクラスを実装してみよう

複数の球で法線レンダリング

これで複数の球を扱えるようになったので、さっそく球を2つ追加して法線表示を行ってみましょう。

#include <iostream>
#include <memory>
#include "vec3.h"
#include "ray.h"
#include "hit.h"
#include "sphere.h"
#include "image.h"
#include "camera.h"
#include "accel.h"


int main() {
    Image img(512, 512);
    Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));
    
    Accel accel;
    accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, 0, 0), 1.0)));  //球1
    accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, -10001, 0), 10000))); //球2

    for(int i = 0; i < img.width; i++) {
        for(int j = 0; j < img.height; j++) {
            double u = (2.0*i - img.width)/img.width;
            double v = (2.0*j - img.height)/img.height;
            Ray ray = cam.getRay(u, v);

            Hit hit;
            if(accel.intersect(ray, hit)) {
                img.setPixel(i, j, (hit.hitNormal + Vec3(1, 1, 1))/2.0); //法線表示
            }
            else {
                img.setPixel(i, j, Vec3(0, 0, 0));
            }
        }
    }

    img.ppm_output();
}
実行結果
実行結果

画像が逆さまになってしまっていますね。でもこれはバグではありません。ピンホールカメラ を使っているので、逆さまの画像が出力されるのです。

しかし逆さまに表示されるのは不便ですね。ピンホールカメラの実装を変更して直しましょう。

Ray getRay(double u, double v) const {
    v = -v; //変更
    Vec3 pinhole = camPos + camForward;
    Vec3 sensorPos = camPos + u*camRight + v*camUp;
    return Ray(camPos, normalize(pinhole - sensorPos));
}
実行結果
実行結果

影の計算

複数の物体を扱えるようになったので、いよいよ影の計算をやってみましょう。

まず前回と同じように、光源の方向l\vec{l}を指定します。

Vec3 lightDir = normalize(Vec3(1, 1, -1))

前回と同じようにレイトレーシングします。

for(int i = 0; i < img.width; i++) {
  for(int j = 0; j < img.height; j++) {
    double u = (2.0*i - img.width)/img.width;
    double v = (2.0*j - img.height)/img.height;
    Ray ray = cam.getRay(u, v);
    Hit hit;
    if(accel.intersect(ray, hit)) {
      //物体に色をつける処理
    }
    else {
      img.setPixel(i, j, Vec3(0, 0, 0));
    }
  }
}

物体上の点p\vec{p}に影がつくということは、点p\vec{p}から光源の方向l\vec{l}を見たときに何らかの物体がそこに存在するので、光が遮られるということです。

影の計算
影の計算

つまり、光源の方向l\vec{l}に物体があるかどうか判定し、物体があれば影をつけるという処理にすればよいのです。

Ray shadowRay = Ray(hit.hitPos, lightDir);
Hit hit2;
//光源との間に物体がある
if(accel.intersect(shadowRay, hit2)) {
  //影をつける
}
else {
  //影をつけない
}

これを実装すると次のようになります。

#include <iostream>
#include <memory>
#include <algorithm>
#include "vec3.h"
#include "ray.h"
#include "hit.h"
#include "sphere.h"
#include "image.h"
#include "camera.h"
#include "accel.h"


int main() {
    Image img(512, 512);
    Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));
    
    Accel accel;
    accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, 0, 0), 1.0)));
    accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, -10001, 0), 10000)));

    Vec3 lightDir = normalize(Vec3(1, 1, -1));

    for(int i = 0; i < img.width; i++) {
        for(int j = 0; j < img.height; j++) {
            double u = (2.0*i - img.width)/img.width;
            double v = (2.0*j - img.height)/img.height;
            Ray ray = cam.getRay(u, v);

            Hit hit;
            //最初に飛ばすRayが当たった場合
            if(accel.intersect(ray, hit)) {
                Ray shadowRay = Ray(hit.hitPos, lightDir); //光源の方向に飛ばすRay
                Hit hit_shadow;
                //もし当たらなかったら色をつける
                if(!accel.intersect(shadowRay, hit_shadow)) {
                    double I = std::max(dot(hit.hitNormal, lightDir), 0.0); //光の強さ
                    img.setPixel(i, j, I*Vec3(1, 1, 1));
                }
                //当たったら影にする
                else {
                    img.setPixel(i, j, Vec3(0, 0, 0));
                }
            }
            else {
                img.setPixel(i, j, Vec3(0, 0, 0));
            }
        }
    }

    img.ppm_output();
}

実行結果
実行結果
ノイズだらけの不気味な画像が出てきました.

なぜこんな結果になってしまったのでしょう。

コンピューターによる浮動小数点数の演算は誤差がつきものという話をしましたね。この結果は実はそこからきています。最初のRayの衝突点を求めるという計算において、わずかな誤差が発生するので、衝突点が床の下側にめり込んでしまっているのです。そのために、そこから光源の方向に伸びるRayが必ず衝突してしまうことになり、影がついて黒くなってしまうのです。

計算誤差
計算誤差

これを防ぐには、衝突点を法線方向にわずかにずらすという処理を加えます。

計算誤差への対処
計算誤差への対処

Ray shadowRay = Ray(hit.hitPos + 0.001*hit.hitNormal, lightDir);
Hit hit2;
//光源との間に物体がある
if(accel.intersect(shadowRay, hit2)) {
  //影をつける
}
else {
  //影をつけない
}

プログラム全体は次のようになります。

#include <iostream>
#include <memory>
#include <algorithm>
#include "vec3.h"
#include "ray.h"
#include "hit.h"
#include "sphere.h"
#include "image.h"
#include "camera.h"
#include "accel.h"


int main() {
    Image img(512, 512);
    Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));
    
    Accel accel;
    accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, 0, 0), 1.0)));
    accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, -10001, 0), 10000)));

    Vec3 lightDir = normalize(Vec3(1, 1, -1));

    for(int i = 0; i < img.width; i++) {
        for(int j = 0; j < img.height; j++) {
            double u = (2.0*i - img.width)/img.width;
            double v = (2.0*j - img.height)/img.height;
            Ray ray = cam.getRay(u, v);

            Hit hit;
            if(accel.intersect(ray, hit)) {
                Ray shadowRay = Ray(hit.hitPos + hit.hitNormal*0.001, lightDir); //法線方向にわずかにずらす
                Hit hit_shadow;
                if(!accel.intersect(shadowRay, hit_shadow)) {
                    double I = std::max(dot(hit.hitNormal, lightDir), 0.0);
                    img.setPixel(i, j, I*Vec3(1, 1, 1));
                }
                else {
                    img.setPixel(i, j, Vec3(0, 0, 0));
                }
            }
            else {
                img.setPixel(i, j, Vec3(0, 0, 0));
            }
        }
    }

    img.ppm_output();
}

実行結果
実行結果
やった!! ちゃんとした影がつきました!!

乱数の生成

ここで急に話の流れが変わりますが、乱数の生成方法について説明しましょう。

C言語の方法

古いC言語の方法では次のようにして乱数を生成していました。

#include <cstdio>
#include <cstdlib>


int main() {
    //シード値
    srand(1);

    for(int i = 0; i < 100; i++) {
        printf("%f\n", (double)rand()/RAND_MAX);
    }
    return 0;
}

実行結果

0.840188
0.394383
0.783099
0.798440
0.911647
0.197551
0.335223
0.768230
...

srand()には シード値(seed) という値をセットします。コンピューターで生成できる乱数はあくまで疑似乱数であり、何らかの漸化式などに従って数列を次々と生成していきます。このとき初期値が必要になるのですが、これがシード値です。

rand()は0からRAND_MAXまでの整数の疑似乱数を返します。プログラムでは、それをRAND_MAXで割ることで、0から1までの疑似乱数にしています。

シード値が同じなら生成される疑似乱数列は全く同じものになります。なので、通常はシード値に時刻などの値をセットします。

#include <cstdio>
#include <cstdlib>
#include <ctime>


int main() {

    //シード値
    srand((unsigned long)time(NULL));

    for(int i = 0; i < 100; i++) {
        printf("%f\n", (double)rand()/RAND_MAX);
    }
    return 0;
}

こうすると実行するたびに表示される疑似乱数列が変わるはずです。

この疑似乱数は 線形合同法 という方法で生成されているのですが、科学技術計算に用いようとする場合には次のような問題点があります。

  • 周期が短い
  • 分布が一様でない

メルセンヌ・ツイスタ

STLにはrandomという乱数を扱うライブラリーが存在します。その中にあるメルセンヌ・ツイスタと呼ばれる疑似乱数生成方法は

  • 周期が長い
  • 分布が一様

という優れた特徴を持つために、科学技術計算において使用される疑似乱数生成方法です。これからレイトレにおいて乱数が必要になった場合にはこのメルセンヌ・ツイスタを使用していきます。

メルセンヌ・ツイスタを用いて乱数を生成するには次のようにします。

#include <cstdio>
#include <random>


std::random_device dev_rnd;
std::mt19937 mt(dev_rnd());
std::uniform_real_distribution<> dist(0, 1);
double rnd() {
    return dist(mt);
}


int main() {
    for(int i = 0; i < 100; i++) {
        printf("%f\n", rnd());
    }
    return 0;
}

少し複雑なので一つ一つ見ていきます。

std::random_device dev_rnd;

std::random_deviceはデバイスの物理的な情報(熱など)から乱数を生成してくれるものです。動作が遅いので、多くの乱数生成に使用することはできません。その代わりに擬似乱数生成器の初期値として使用します。

std::mt19937 mt(dev_rnd())

std::mt19937はメルセンヌ・ツイスタによる擬似乱数生成器です。初期値をデバイス乱数で初期化しています。

std::uniform_real_distribution<> dist(0, 1);

std::uniform_real_distribution<>は疑似乱数生成器による乱数を指定された範囲の一様分布に従うように変換してくれるものです。この場合は0から1までの一様分布になるように変換されます。

double rnd() {
  return dist(mt);
}

std::uniform_real_distribution<>の引数に疑似乱数生成器を入れてあげれば、0から1までの一様分布に従う乱数が得られます。

スーパーサンプリング(SSAA)

先ほどのレンダリング画像をもう一度見てみましょう。

先ほどのレンダリング画像
先ほどのレンダリング画像
物体の境界、影の境界がギザギザしていますね。これは本来連続的な曲線である物体の境界などを、離散的な画素で表現しているために現れるもので、エイリアシング(Aliasing) と呼ばれています。

これを防ぐ方法が スーパーサンプリング(SSAA) と呼ばれるものです。

SSAA
SSAA

一つの画素内で始点をわずかにずらして複数のRayを飛ばし、平均した色をその画素に書き込みます。こうすることで、物体の境界線などがやんわりとした感じになり、美しく見えるようになります。

メルセンヌ・ツイスタを使ってこれを実装すると次のようになります。

#include <iostream>
#include <memory>
#include <algorithm>
#include <random>
#include "vec3.h"
#include "ray.h"
#include "hit.h"
#include "sphere.h"
#include "image.h"
#include "camera.h"
#include "accel.h"


std::random_device dev_rnd;
std::mt19937 mt(dev_rnd());
std::uniform_real_distribution<> dist(0, 1);
double rnd() {
    return dist(mt);
}


int main() {
    Image img(512, 512);
    Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));
    
    Accel accel;
    accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, 0, 0), 1.0)));
    accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, -10001, 0), 10000)));

    Vec3 lightDir = normalize(Vec3(1, 1, -1));

    for(int k = 0; k < 100; k++) {
        for(int i = 0; i < img.width; i++) {
            for(int j = 0; j < img.height; j++) {
                //ピクセル内の位置をランダムにずらす
                double u = (2.0*(i + rnd()) - img.width)/img.width;
                double v = (2.0*(j + rnd()) - img.height)/img.height;
                Ray ray = cam.getRay(u, v);

                Vec3 color;
                Hit hit;
                if(accel.intersect(ray, hit)) {
                    Ray shadowRay = Ray(hit.hitPos + hit.hitNormal*0.001, lightDir);
                    Hit hit_shadow;
                    if(!accel.intersect(shadowRay, hit_shadow)) {
                        double I = std::max(dot(hit.hitNormal, lightDir), 0.0);
                        color = I*Vec3(1, 1, 1);
                    }
                    else {
                        color = Vec3(0, 0, 0);
                    }
                }
                else {
                    color = Vec3(0, 0, 0);
                }

                img.setPixel(i, j, img.getPixel(i, j) + 1/100.0*color); //平均をとる
            }
        }
    }

    img.ppm_output();
}

実行結果
実行結果
先ほどよりだいぶ綺麗になりましたね!!

並列化計算

SSAAもやるようになってくると、だいぶ計算に時間がかかるようになってきます。もう少し早く計算できるようになると嬉しいですね。

現在のCPUにはたいてい複数のコアが搭載されているのですが、今まで書いてきたプログラムは一つのコアのみを用いて計算を行っています。

シングルコア
シングルコア

早く計算できるようにするために、CPUのコアを全て使って計算を行うことを考えてみましょう。

マルチコア
マルチコア

このようにfor文を幾つかの領域に分けて、それぞれの領域を別々のCPUコアで計算することができるようになります。

並列計算
並列計算

openMP

openMPは並列化計算を簡単に実現できるようにしてくれるプラットフォームです。

openMPのインストール

https://qiita.com/kaityo256/items/ae9329dae24ea8828ae0

openMPの使い方

int main() {
  for(int i = 0; i < 100; i++) {
    c[i] = a[i] + b[i];
  }
}

並列化

#include <omp.h> //これが必要

int main() {
  #pragma omp parallel for
  for(int i = 0; i < 100; i++) {
    c[i] = a[i] + b[i];
  }
}

#pragma omp parallel forを並列化したいfor文の直前につけるだけです。

さらにコンパイルするときに次のようなオプションを追加してあげます。

g++ -fopenmp main.cpp

これでCPUのコアを全て使って並列化計算できるようになります。

並列化計算の問題点

並列化計算が正しく行われるためには、それぞれの計算タスクが独立 している必要があります。独立していない場合、計算が行われる順序によって計算結果が変わってきてしまいます。

しかしそれを防ごうとして計算順序を指定すると、その計算が終わるまで他の計算を進めることができなくなるので、結果的に並列化の恩恵を受けられなくなります。

つまり、何でもかんでも並列化すればいいというわけではなく、計算タスクが独立している部分だけに使用するべきです。

幸い、レイトレーシングは画素ごとに色を求める計算は互いに独立しているので、この部分に対して並列化を適用することができます。

並列化計算によるレイトレ

#include <iostream>
#include <memory>
#include <algorithm>
#include <random>
#include <omp.h>
#include "vec3.h"
#include "ray.h"
#include "hit.h"
#include "sphere.h"
#include "image.h"
#include "camera.h"
#include "accel.h"


std::random_device dev_rnd;
std::mt19937 mt(dev_rnd());
std::uniform_real_distribution<> dist(0, 1);
double rnd() {
    return dist(mt);
}


int main() {
    Image img(512, 512);
    Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));
    
    Accel accel;
    accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, 0, 0), 1.0)));
    accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, -10001, 0), 10000)));

    Vec3 lightDir = normalize(Vec3(1, 1, -1));

#pragma omp parallel for
    for(int k = 0; k < 100; k++) {
        for(int i = 0; i < img.width; i++) {
            for(int j = 0; j < img.height; j++) {
                //ピクセル内の位置をランダムにずらす
                double u = (2.0*(i + rnd()) - img.width)/img.width;
                double v = (2.0*(j + rnd()) - img.height)/img.height;
                Ray ray = cam.getRay(u, v);

                Vec3 color;
                Hit hit;
                if(accel.intersect(ray, hit)) {
                    Ray shadowRay = Ray(hit.hitPos + hit.hitNormal*0.001, lightDir);
                    Hit hit_shadow;
                    if(!accel.intersect(shadowRay, hit_shadow)) {
                        double I = std::max(dot(hit.hitNormal, lightDir), 0.0);
                        color = I*Vec3(1, 1, 1);
                    }
                    else {
                        color = Vec3(0, 0, 0);
                    }
                }
                else {
                    color = Vec3(0, 0, 0);
                }

                img.setPixel(i, j, img.getPixel(i, j) + 1/100.0*color); //平均をとる
            }
        }
    }

    img.ppm_output();
}

並列計算が正しく動いていれば、計算時間が先ほどより短くなっているはずです。

アンビエントオクルージョン(AO)

ある点から空がどれくらいの割合で見えるかを表すのが アンビエントオクルージョン(AO) です。物体が隣接している場所は空が見える範囲が狭くなるので割合が小さくなり、逆に何もないところでは0.5に近づきます。この値を可視化してみましょう。

ランダムな方向の生成

ランダムに単位球上の点への方向を生成する関数Vec3 randomDir()を作ってみましょう。

天球全体へのランダムな方向を作るには、単位球上の点をランダムに生成して、原点からその点への方向を代わりに使えば良いです。

天球全体へのランダムな方向生成
天球全体へのランダムな方向生成

そこで、単位球上の点をランダムに生成することを考えます。球面座標系(ϕ,θ)(\phi, \theta)ϕ\phiθ\thetaをランダムに生成すると、単位球上の点(x,y,z)(x, y, z)

x=cosϕsinθ   y=cosθ   z=sinϕsinθx = \cos{\phi}\sin{\theta} ~~~ y = \cos{\theta} ~~~ z = \sin{\phi}\sin{\theta}

となります。これをそのまま方向ベクトルとして使えばよいです。(実は一様サンプリングにはなっていない)

Vec3 randomDir() {
    double theta = M_PI*rnd();
    double phi = 2*M_PI*rnd();
    
    double x = std::cos(phi)*std::sin(theta);
    double y = std::cos(theta);
    double z = std::sin(phi)*std::sin(theta);
    return Vec3(x, y, z);
}

AOの計算

あとはshadowRayの方向をrandomDir()で生成し、結果の平均を取るようにすればよいです。

#pragma omp parallel for
for(int k = 0; k < 100; k++) {
    for(int i = 0; i < img.width; i++) {
        for(int j = 0; j < img.height; j++) {
            //ピクセル内の位置をランダムにずらす
            double u = (2.0*(i + rnd()) - img.width)/img.width;
            double v = (2.0*(j + rnd()) - img.height)/img.height;
            Ray ray = cam.getRay(u, v);

            Vec3 color;
            Hit hit;
            if(accel.intersect(ray, hit)) {
                Vec3 lightDir = randomDir(); //ランダムな方向を生成
                Ray shadowRay = Ray(hit.hitPos + hit.hitNormal*0.001, lightDir);
                Hit hit_shadow;
                if(!accel.intersect(shadowRay, hit_shadow)) {
                    color = Vec3(1, 1, 1);
                }
                else {
                    color = Vec3(0, 0, 0);
                }
            }
            else {
                color = Vec3(1, 1, 1);
            }

            img.setPixel(i, j, img.getPixel(i, j) + 1/100.0*color); //平均をとる
        }
    }
}

実行結果
実行結果
影がとてもリアルに見えますね!!

AOの計算でやっていることを数式で表す

天球全体をΩ\Omega、点xxから方向ω\vec{\omega}が見えるということをV(x,ω)V(x, \vec{\omega})で表すと、点xxにおけるAOの計算は次のように表すことができます。

AO=14πΩV(x,ω)dωAO = \frac{1}{4\pi}\int_{\Omega} V(x, \vec{\omega}) d\vec{\omega}

V(x,ω)V(x, \vec{\omega})xxからω\vec{\omega}の方向に何もないときに11、それ以外のときは00になる関数です。方向ω\vec{\omega}のことを 立体角(Solid Angle) といいます。全球の立体角は4π4\piになるので、積分を4π4\piで割っています。

上の計算では、ランダムな方向列{ω1,ω2,,ωn}\{ \vec{\omega_1}, \vec{\omega_2}, \dots, \vec{\omega_n} \}を生成し、

AO=1ni=1nV(x,ωi)AO = \frac{1}{n}\sum_{i = 1}^n V(x, \vec{\omega_i})

とすることで、上の積分を近似していました。一般的には上の積分と、この近似式の値はnnを大きくしても一致することはありません。

上の積分を正しく近似して計算した結果は以下のようになります。

正しい計算結果
正しい計算結果
先ほどの近似計算による結果とはかなり異なります。

では、上の積分に一致するように近似式を作るにはどうすればよいのでしょうか?この続きは物理ベースレンダリング理論のところでやります。

練習問題

問1

絶対屈折率n1n_1から絶対屈折率n2n_2の媒質に入射するレイd\vec{d}がある。媒質境界面の法線をn\vec{n}としたときに、レイが屈折した方向r\vec{r}を求めよ。なおn1<n2n_1 < n_2とする。

屈折方向
屈折方向

補足

スネルの法則により、入射角と屈折角には n1sinθi=n2sinθtn_1\sin{\theta_i} = n_2\sin{\theta_t} の関係がある。

問2

講義でランダムな方向を生成する関数randomDir()を作ったが、実はこの関数は単位球面上の点の一様サンプリングにはなっていない。例えば、北極、南極点に近いところでは点が密集する傾向にある。単位球面上での一様サンプリングを行うにはどうすればよいだろうか?計算方法を考えてみよ。

問3

円周率を乱数を用いて評価してみよう。 [0,2]×[0,2][0, 2]\times[0, 2]の正方形内の点をnn個ランダムに生成する。これらの点が中心座標(1,1)(1, 1)で半径11の円に入った回数をaaとする。すると

π4an\pi \approx 4\frac{a}{n}

で推定することができる。この方法で円周率を求めるプログラムを作成せよ。