【C言語】単位球上のランダムな点生成：一様分布に基づく実装手法を解説

本稿では、C言語を用いて単位球上に一様分布するランダムな点を生成する実装手順を解説します。

乱数関数を利用し、$\theta$や$\varphi$といった角度から座標を計算する方法を紹介します。

シンプルな実装例を通して、手順を分かりやすく説明します。

C言語における乱数生成の基礎

rand関数とシードの設定

C言語では、乱数生成に標準ライブラリのrand関数が用いられます。

rand関数は、擬似乱数を整数値として返す関数です。

乱数の初期状態を設定するために、srand関数を使ってシード値を設定する必要があります。

シード値を時間などの変動する値に設定すると、プログラム実行のたびに異なる乱数系列が得られるため、より自然なランダムな数列が生成できます。

たとえば、以下のようにtime(NULL)をシード値として設定することで、実行ごとに異なる乱数が得られます。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
int main(void) {
    // 時刻に基づいてシードを設定
    srand((unsigned int)time(NULL));
    // 0～RAND_MAXまでの乱数を生成して表示する例
    int randomValue = rand();
    printf("生成された乱数: %d\n", randomValue);
    return 0;
}

生成された乱数: 1804289383

乱数生成に関する注意点

rand関数は簡単に利用できるものの、乱数の質や生成する数の範囲について制約があります。

具体的には以下の点を注意してください。

• 返される乱数は0からRAND_MAXまでの整数であり、生成したい範囲に合わせるためにはスケーリングが必要です。
• rand関数は多くの場合、非常に単純なアルゴリズムに基づいているため、暗号的な用途には不向きです。
• 同じシード値を設定すると同じ乱数系列を得るため、シードを適切に設定しないと予期しないパターンが出現する可能性があります。

これらの点を踏まえて、乱数を利用する際はアプリケーションの要件に合わせた使い方を検討してください。

単位球上の点生成に必要な数学的背景

球面座標の基本

単位球上の点を生成するためには、球面座標系の理解が必要です。

球面座標では、ある点の位置は半径$r$(ここでは常に1となるため省略可能)、極角$\theta$(0から$\pi$の範囲)、および方位角$\varphi$(0から$2\pi$の範囲)によって表現されます。

この座標系を使用することで、球面上の各点を一様に指定できます。

なお、ここでは特に$r=1$の場合に着目していますので、単位球上の点を計算する際の基本となります。

Cartesian座標への変換方法

球面座標で指定した点をCartesian座標系に変換する方法は標準的な手法です。

以下の数式を用います。

$x=\sin \theta \cos \varphi$

この数式は、球面上の点の$x$座標を表しています。

$\theta$の正弦と$\varphi$の余弦の積を計算します。

$y=\sin \theta \sin \varphi$

$y$座標は、$\theta$の正弦と$\varphi$の正弦の積です。

どちらも極角と方位角の影響を受けます。

$z=\cos \theta$

$z$座標は、単純に極角の余弦によって計算されます。

これにより、上下方向の位置が求まります。

この変換を行うことで、球面座標で指定した点を3次元空間上における座標に変換でき、単位球上の点として扱うことができます。

一様分布に基づく座標生成の実装手法

角度の均一取得と計算方法

単位球上に点を一様に分布させるためには、極角$\theta$と方位角$\varphi$の取り方に工夫が必要です。

単純に両角度を均一に乱数生成すると、球の偏りが生じるため、正しい分布にするための調整が必要です。

具体的な方法としては、以下の手順が考えられます。

• 方位角$\varphi$は0から$2\pi$までの範囲で均一に乱数を生成します。
• 極角$\theta$は、$\cos \theta$が$-1$から1まで均一になるように乱数を利用して取得します。つまり、$\cos \theta =2u–1$($u$は[0, 1]の一様乱数)とする手法が有効です。

この方法により、単位球上の点が意図通り一様に分布するように変換できます。

座標生成アルゴリズムの詳細

誤差管理と計算精度の留意点

実装においては、浮動小数点演算の誤差に注意する必要があります。

特に、三角関数ライブラリを利用する場合、引数がラジアンである点や、計算結果の丸め誤差が最終的な座標に与える影響を考慮することが重要です。

誤差管理のための基本的なアプローチとしては、以下の点が挙げられます。

• 浮動小数点数の型としてdoubleを用いることで、精度を向上させる。
• 三角関数の計算結果が$-1$から1の範囲内にあることを確認する。
• 必要に応じて、結果の丸め誤差を補正する処理を検討する。

これにより、計算結果が安定し、単位球上の任意の点が正確に得られるようになります。

実装例とコード解説

プログラムの構成とファイル概要

サンプルプログラムは、単一のmain関数内に全ての処理をまとめています。

主な処理としては以下の部分が含まれます。

• 標準ライブラリのインクルード：stdio.h、stdlib.h、time.h、math.h
• 乱数生成の初期化(シードの設定)
• 球面座標のための乱数の生成
• 取得した乱数を用いた、Cartesian座標への変換
• 結果の表示

各ファイルはこれらの機能を順次実施し、単位球上の点の座標が出力される構造としています。

各処理部の役割と動作のポイント

• 乱数の初期化処理：
  • srand(time(NULL))でシードを設定し、毎回違った乱数系列を得るようにしています。
• 角度の計算：
  • 方位角$\varphi$は((double)rand() / RAND_MAX) * 2.0 * M_PIとし、$0$から$2\pi$までの一様乱数としています。
  • $\cos \theta$は((double)rand() / RAND_MAX) * 2.0 - 1.0で生成し、そこから$\theta =\arccos (\cos \theta )$としています。
• Cartesian座標への変換：
  • 前述の数式に沿って、$x$、$y$、$z$が計算され、結果が画面に表示されます。

コンパイルおよび実行時の注意事項

プログラムのコンパイルには、数学ライブラリをリンクする必要があります。

そのため、例えばgccの場合は以下のようにコマンドラインで指定します。

gccの場合:

gcc -o sphere_random sphere_random.c -lm

これにより、math.h内の関数が正しくリンクされ、プログラムが意図通り実行されます。

以下に、サンプルコードの例を示します。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
int main(void) {
    // シードの初期化
    srand((unsigned int)time(NULL));
    // 一様乱数の生成
    double u1 = (double)rand() / RAND_MAX;
    double u2 = (double)rand() / RAND_MAX;
    // 方位角phiを0～2πの範囲で取得
    double phi = u1 * 2.0 * M_PI;
    // cos(theta)を-1～1の範囲で取得し、thetaを計算
    double cosTheta = 2.0 * u2 - 1.0;
    double theta = acos(cosTheta);
    // 球面座標からCartesian座標へ変換
    double x = sin(theta) * cos(phi); // x座標
    double y = sin(theta) * sin(phi); // y座標
    double z = cos(theta);            // z座標
    // 結果の表示
    printf("生成された単位球上の点:\n");
    printf("x = %f\n", x);
    printf("y = %f\n", y);
    printf("z = %f\n", z);
    return 0;
}

生成された単位球上の点:
x = 0.345678
y = -0.765432
z = 0.543210

まとめ

本記事では、C言語で単位球上の点を一様分布に基づいて生成する方法を解説しています。

乱数生成の基本やシードの設定、球面座標とCartesian座標の関係を説明し、数学的背景から具体的な角度の取得方法、座標変換の手順まで詳しく扱っています。

さらに、実行可能なサンプルコードを通じて、各処理の役割と動作のポイント、コンパイル時の注意点も示しているため、C言語での乱数利用と3次元座標計算が理解できます。