エフアンダーバー

個人でのゲーム開発

乱数生成の備忘録2

他人が書いた乱数生成のコードを読んでいたところ、 前回自分が書いたものと違う実装をしていたので改めていろいろと調査。

www.f-sp.com

問題のコード

異なる実装となっていたのは整数の乱数から[0,1)の範囲の実数の乱数を生成するコード。

前回自分が書いたコードでは次のように定数を掛けて変換していた。

static Double ToDouble(UInt64 x)
{
    return (x * ((1.0 - 1.0 / (1ul << 52)) / UInt64.MaxValue));
}

今回見つけたコードは次のような実装となっていた*1

static Double ToDouble(UInt64 x)
{
    const UInt64 ExponentBits = 0x3ff0_0000_0000_0000;
    const UInt64 MantissaMask = 0x000f_ffff_ffff_ffff;
    UInt64 bits = (x & MantissaMask) | ExponentBits;
    Double v = BitConverter.Int64BitsToDouble(unchecked((Int64)bits));
    return (v - 1.0);
}

double型(倍精度浮動小数点数)はIEEE754規格に従っているので、 1bitの符号部と11bitの指数部および52bitの仮数部から成る。

この実装では符号部を正(0)とし、指数部が0(3ff)となるようにして、残りの仮数部を乱数の値で埋めている。 このようにすると、数値は単純に仮数部の先頭に1を加えた1.XXX...XXXで表される[1,2)の範囲の値となる。 最後にその値から1を引けば値は[0,1)の範囲に収まることになる。

メジャーなプロジェクトとしてV8のコードを見てみたところ、こちらも後者で実装されていた。

何が違うか

前者のコードでは一様分布な整数の乱数の範囲を[0,1)に調整した後で値が丸められる。 異なる実数であっても、丸めによって同じ値となってしまうことがあるため、結果の浮動小数点数自体の分布は一様ではなくなる。 浮動小数点数が表現できる値には偏りがあり、0に近いほど高い精度で表現できるため、 比較的精度の低い1付近では0付近に比べて丸めによる同じ値の発生が起きやすい。 ただし、表現している実数自体は一様に分布している。

一方、後者のコードは計算機イプシロンの間隔で並ぶ[0,1)の範囲の浮動小数点数の部分集合のみを利用して、それらに一様の分布を与えている。 つまり、[0,1)の範囲で最も精度の悪い部分に合わせることにより、浮動小数点数(の部分集合)自体の分布をも一様としている。 欠点は[0,1)の範囲を表現する浮動小数点数62bit分の内の52bit分(1/1024)しか使用できないこと。 しかし、削られた数値からの距離が計算機イプシロンの半分以下である値が必ず含まれるので気にするほどではないかもしれない。

どちらがいいか

いくらか調べてはみたもののどちらかが決定的にダメだと思えるような例は見当たらなかった。 ただし、後者を強く支持する人がいる一方で、前者はこれも悪いとは言えないよねくらいの感じだったので後者の方が無難かもしれない。

しかし、後者は実装する言語によっては少し面倒くさい。 例えば、C#ではBitConverter.Int64BitsToDoubleに対応するfloat(単精度浮動小数点数)用のメソッドが見当たらない。 そのため次のように少し回りくどいことをしなければならない。

  • byte配列を介してBitConverterで変換
    • ヒープメモリの確保が必要
    • オブジェクトを使いまわすならスレッドセーフでなくなる(あるいはThreadStaticAttributeが必要)
  • unsafeでreinterpret_cast的な変換
  • StructLayoutAttributeで共用体的な変換

・・・とか考えていたところ、最初に見つけたコード(JavaScript)が次のようなコードに更新されていた(実際には定数部分は最適化されている)。

function toDouble(upper, lower) {
    return upper * Math.pow(2, -32) + (lower >> 12) * Math.pow(2, -52);
}

C#で書くならこんな感じ。

static Double ToDouble(UInt64 x)
{
    return (x >> 12) * (1.0 / (1ul << 52));
}

これならfloat(Single)でも同じように書ける。

static Single ToSingle(UInt32 x)
{
    return (x >> 9) * (1.0f / (1u << 23));
}

乗算が含まれるとはいえ、回りくどい変換をするよりはこちらの方が速度も良さそうな気がする。

参考

おわりに

調べてみて改めて思ったのは乱数とか浮動小数点数とかすごくメンドクサイな、と。 正しいか正しくないかの検証が楽ならいいのだけど、どうなれば正しいのかがよくわからないのでいまいち結果に自信を持てない。

書いている途中で気づいたのだけど、前回の記事でdouble(64bit)の乱数を生成するのにuint(32bit)の乱数を元にしているという。 分布がどうのこうの言うより前にまともなコードを書かねば・・・。

*1:下記は同様の実装となるようにC#で書き換えたもの

乱数生成の備忘録

以前書いた疑似乱数生成コードの修正作業をしたので、その際に調べたことや考えたことのメモ。

基本的に自分用なので雑。

動機と方針

  • System.Randomのアルゴリズムは古めで偏りがある。
    • Knuthの減算法。
    • 再現性の問題があるので今後も変わることはなさそう。
  • UnityEngine.Randomなどゲームエンジンや外部ライブラリに依存したくない。
  • System.Randomは継承しない。
    • .NET的には推奨らしいけど無駄なフィールドのメモリが必要なため。

0~1の範囲の乱数

1を含めるか否か

各APIの仕様は、

  • System.Randomは1を含めない。
  • JavaScriptのMath.randomも1を含めない。
  • UnityEngine.Randomは1を含める。

利用側として、

  • 1を含めて困るケースは考えられる。
    • 適当な整数を掛けてから整数にキャストした場合に滅多に現れないが現れうる数が発生。
      • 整数の乱数は別にとれるのでキャストの必要があるかどうかは微妙。
  • 1を含めないと困るケースは思いつかない。

値域は[0,1)の方が無難に思える。

最大値は何か

1を含めないとなると最大値として何をとるべきか。

リファレンスによるとSystem.RandomのNextDoubleでは最大値として0.99999999999999978を返すらしい。 いろいろと試した結果、どうやら計算機イプシロンをeとして(1-e)の値っぽい。

Console.WriteLine("{0:E16}", 1.0 - Math.Pow(2, -52));
// 9.9999999999999978E-001

乱数は次のメソッドでUInt32の範囲の値を返すようにしたので、

public abstract UInt32 Sample();

SingleとDoubleの乱数を返すメソッドはそれぞれ次のように実装。

public Single NextSingle()
{
    return (Sample() * ((1.0f - 1.0f / (1u << 23)) / UInt32.MaxValue));
}

public Double NextDouble()
{
    return (Sample() * ((1.0 - 1.0 / (1ul << 52)) / UInt32.MaxValue));
}

誤差やら精度やらの問題があるのでいまいち自信ないけどとりあえずはこれで。

最大値は妥当か

MDNのrandomの説明を読むと気になることが書いてある。

JavaScript における数値は、IEEE 754 浮動小数点での最近接の偶数への丸め (round-to-nearest-even) の振る舞いをするので、 以下の関数についての値域 (ただし Math.random() それ自身についての値域は除く) が厳密ではない、ということに注意してください。 非常に大きい境界値(2^53 以上)を選ぶと、極めて稀なケースにおいては、通常なら除外される上限値を算出してしまうことがあり得ます。

randomもSystem.RandomのNextDoubleと同じような実装をしていると思われるので、上述の実装でも同じことが起こるのかもしれない。

しかし、非常に大きい数を掛けることはほぼないと思われるのでとりあえず無視しておくことにする。

Int32からUInt64への変換

以前書いたコードを読んでいたところ気になるコードを発見。

UInt64 y = unchecked((UInt32)x); // Int32 x;

UInt64へのキャストの間違いかなと思いつつ、一応テスト。

Console.WriteLine(unchecked((UInt64)(-1)));         // 18446744073709551615
Console.WriteLine(unchecked((UInt64)(UInt32)(-1))); // 4294967295

どうもInt32から直接UInt64へとキャストすると、ゼロ拡張ではなく符号拡張になる模様。 UInt32経由ではなくInt64経由でキャストした挙動になるみたい。

ちなみに上のコード、(-1)の括弧を外すとマイナス符号が単項演算子なのか二項演算子なのか判断しかねるらしく、「ulongは型だろ」と怒られる(知ってた)。

さらにはその際の日本語のエラーメッセージは

error CS0119: ‘ulong’ は 種類 です。これは特定のコンテンツでは無効になります。

原文は

error CS0119: ‘ulong’ is a type, which is not valid in the given context.
(‘ulong’ は 型 です。これはこの文脈では有効ではありません。)

全然違うじゃん…

おわりに

まともな記事として書けるほど内容もなかったので、とりあえずメモ程度にまとめ。

メモなので、いつものですます調も変かと思って封印してみたけど、使えないとなるとそれもそれで書きづらかったり。 まあ、文章の練習だと思って頑張ろう。

なんか記事によって口調違うと情緒不安定な人みたいだな…

【Unity】 Rigidbodyの回転方法

Rigidbody(剛体)の移動方法については以前書きましたが、回転方法については書いていなかったので今更ながらまとめてみます。

まあ、移動方法とそんなに変わらないんですが・・・。

www.f-sp.com

準備

説明の簡略化のため、本記事内では下記の記号は特別な記述がない限り対応する物理量を表すものとします。

記号 物理量
t 時刻
\omega 角速度ベクトル
\alpha 角加速度ベクトル
I 慣性モーメントテンソル
N トルク

物理量の意味がわからない場合には、記事末尾の付録を確認してください。

Rigidbodyの回転

Rigidbodyの基本的な回転方法には次の四つがあります。

  • rotation
  • angularVelocity
  • AddTorque
  • MoveRotation

以降のサンプルコードでは度々rigidbodyという変数が登場しますが、 これはMonoBehaviourのrigidbodyプロパティではありません(非推奨のプロパティ)。 事前に次のようなコードでrigidbodyを取得してください。

Rigidbody rigidbody = GetComponent<Rigidbody>();

rotation

rotationはRigidbodyの姿勢を変更します。

rigidbody.rotation = Quaternion.AngleAxis(180.0f, Vector3.up);

Rigidbodyをアタッチしたゲームオブジェクトの姿勢を変更したい場合には、 Transform.rotationではなく、Rigidbody.rotationを変更します。

この方法による姿勢の変更は変化ではなく上書きです。 要するに、元々どんな姿勢だったかは一切考慮されません。 そのため、初期姿勢の設定には適していますが、特定の方向に向き直るのには不適です。 変化を表したいときにはMoveRotationを利用します。

angularVelocity

angularVelocityはRigidbodyの角速度ベクトルを変更します。

rigidbody.angularVelocity = Vector3.up * Mathf.PI;

元々の角速度ベクトルを上書きしていることに注意が必要です。

AddTorque

AddTorqueはトルク(力のモーメント)を加えます。角加速度ベクトルを設定すると考えてもいいかもしれません。

rigidbody.AddTorque(Vector3.up * Mathf.PI);

AddTorqueは第二引数にForceModeというオプションがあり、それによってトルクの加え方が変わります(省略した場合はForceMode.Forceです)。

ForceMode 質量 タイプ
Force 考慮 継続的
Impulse 考慮 瞬間的
Acceleration 無視 継続的
VelocityChange 無視 瞬間的

このForceModeは移動の際のAddForceのオプションと同じです。 トルクは力のモーメントであり、トルクを加えるというのは偶力(平行で大きさの等しい逆向きの二つの力)を加えるようなものです。 トルクという概念で扱ってはいますが、実際には力を加えているのでAddForceと同じオプションが利用できるのです。

ForceMode.Force

Forceは物理における「力」です。 何かを押し続けるというような継続的な力を表します。 『プレイヤーがキーを押している間』や『ある範囲にいる間』というように特定の期間中、トルクを与え続ける場合に使います。

AddTorqueの第一引数  N_{Force} I\alpha あるいは  \frac{I\Delta\omega}{\Delta t} を表す値と解釈されます。

I_p を主慣性モーメント、R を慣性主軸からローカル座標軸への回転とすると、

 \displaystyle
    \Delta\omega = (I^{-1} N_{Force})\Delta t = (RI_p^{-1}R^{-1} N_{Force})\Delta t

なので次の二行のコードは同じような意味になります *1

rigidbody.AddTorque(Vector3.up * Mathf.PI, ForceMode.Force);

var I = rigidbody.inertiaTensor;
var R = rigidbody.rotation * rigidbody.inertiaTensorRotation;
rigidbody.angularVelocity += R * Vector3.Scale(Reciprocal(I), Quaternion.Inverse(R) * (Vector3.up * Mathf.PI)) * Time.fixedDeltaTime;

Reciprocalは次のように定義されているものとします。

static Vector3 Reciprocal(Vector3 v)
{
    return new Vector3(1.0f / v.x, 1.0f / v.y, 1.0f / v.z);
}

ForceMode.Impulse

Impulseは物理における「力積」です。 釘を打つときのような瞬間的な力を表します。 『プレイヤーがキーを押した瞬間』や『物体同士が接触した瞬間』というように何かのイベントの発生に対して、一瞬だけトルクを加える場合に使います。

AddTorqueの第一引数  N_{Impulse} I\Delta\omega を表す値と解釈されます。

 \displaystyle
    \Delta\omega = I^{-1} N_{Impulse} = RI_p^{-1}R^{-1} N_{Impulse}

なので次の二行のコードは同じような意味になります。

rigidbody.AddTorque(Vector3.up * Mathf.PI, ForceMode.Impulse);

var I = rigidbody.inertiaTensor;
var R = rigidbody.rotation * rigidbody.inertiaTensorRotation;
rigidbody.angularVelocity += R * Vector3.Scale(Reciprocal(I), Quaternion.Inverse(R) * (Vector3.up * Mathf.PI));

ForceMode.Acceleration

Accelerationは物理における「加速度」です。  \mbox{力} = \mbox{質量} \times \mbox{加速度} なのでForceMode.Forceの質量(及び慣性モーメント)を無視したバージョンになります。

AddTorqueの第一引数  N_{Acceleration} \alpha あるいは  \frac{\Delta\omega}{\Delta t} を表す値と解釈されます。

 \displaystyle
    \Delta\omega = N_{Acceleration}\Delta t

なので次の二行のコードは同じような意味になります。

rigidbody.AddTorque(Vector3.up * Mathf.PI, ForceMode.Acceleration);

rigidbody.angularVelocity += (Vector3.up * Mathf.PI) * Time.fixedDeltaTime;

ForceMode.VelocityChange

VelocityChangeは物理における「速度変化」です。  \mbox{力積} = \mbox{質量} \times \mbox{速度変化} なのでForceMode.Impulseの質量(及び慣性モーメント)を無視したバージョンになります。

AddTorqueの第一引数  N_{VelocityChange} \Delta\omega を表す値と解釈されます。

 \displaystyle
    \Delta\omega = N_{VelocityChange}

なので次の二行のコードは同じような意味になります。

rigidbody.AddTorque(Vector3.up * Mathf.PI, ForceMode.VelocityChange);

rigidbody.angularVelocity += (Vector3.up * Mathf.PI);

MoveRotation

MoveRotationはほぼrotationと同じで、Rigidbodyの姿勢を変更します。

rigidbody.MoveRotation(Quaternion.AngleAxis(180.0f, Vector3.up));

ただし、isKinematicの値によっては微妙に異なる挙動をします (以下、リファレンスに詳細が書かれていないため、リファレンスの内容と実際の挙動からの推測です)。

まず、isKinematicがfalseのときにはrotationと全く同じ挙動となります。 つまり、次の二行のコードは同じ意味になります。

rigidbody.rotation = Quaternion.AngleAxis(180.0f, Vector3.up);

rigidbody.MoveRotation(Quaternion.AngleAxis(180.0f, Vector3.up));

次に、isKinematicがtrueのときにはrotationの上書きとは異なり、変化を表します。 つまり、元の姿勢から指定した姿勢へと瞬時に回転したとみなされます。 これはRigidbodyの設定によってはrotationと挙動が異なります。 例えば、interpolationを設定していた場合には回転の中間状態が補間されてレンダリングされる可能性があります。

実際、処理がどのように違うかというと、 rotationを変更した場合にはそのままrotationの値を変更しているだけなのですが、 MoveRotationを呼び出した場合には回転量と回転時間から角速度を逆算して、 その角速度で回転するように設定しているようです。 isKinematicがtrueの状態なので次の更新処理時には角速度は再度0に戻ります。

その他

前述の四つの回転方法はどれも直接回転パラメータを指定していましたが、 回転は本来、重心から離れた点に横向きの力が加わることで発生します。 摩擦でRigidbodyが意図しない回転をしてしまうのはよくあるバグだと思いますが、 これも重心から離れた表面に接線方向の力が働いたことが原因です。

現実的な回転をさせたい場合には、 移動と回転それぞれを指定するのではなく、 AddForceAtPositionによって位置を指定して力を与えるのが効果的です。 AddForceAtPositionの使い方は位置を指定する以外はAddForceと同じなので過去記事を参照してください。

付録:剛体と回転

Rigidbodyは日本語では「剛体」と呼ばれる物理学上のモデルです。

剛体は力を加えても一切変形しない物体で、質点がたくさん集まって形を成したものとして扱われます。 質点の持つ「位置」や「質量」といった概念に、 新たに「形状*2」という概念が加わったものくらいに考えるといいと思います。

大きさや形のない質点では物体の位置しか考えませんでしたが、 剛体には形状があるため、現在どんな姿勢をとっているかについて考える必要があります。 剛体の位置に関する運動を「並進運動」、姿勢に関する運動を「回転運動」と言います。

剛体の回転運動では、並進運動に用いる様々な物理量に対応した物理量を考えます。 わかりやすいところで言うと「速度」に対応する「角速度」なんかがあります。

並進運動 表記 回転運動 表記
位置  r 角度  \theta
速度  v = \frac{dr}{dt} 角速度  \omega = \frac{d\theta}{dt}
加速度  a = \frac{dv}{dt} 角加速度  \alpha = \frac{d\omega}{dt}
(慣性)質量  m 慣性モーメント  I
運動量  p = mv 角運動量  L = r \times p
 F トルク  N = r \times F

そして、素晴らしいことにそれらに対して並進運動とよく似た関係が成り立ちます。 例えば、ニュートンの運動方程式 F = maのように N = I\alphaが成り立ちます。

角速度・角加速度

並進運動における速度・加速度に対応する値として、 回転運動には角速度・角加速度があります。 これらは通常ベクトルで表されるため、それぞれ角速度ベクトル・角加速度ベクトルと言ったりもします。

回転運動はある一瞬で見れば常にひとつの平面上で行われます。 そのため、どの平面でどれくらいの角度の回転をさせるかの情報があれば三次元の回転を表現できます。 ベクトルはこれらの情報を保持するのに便利で、ベクトルの向きでそれに垂直な平面を、ベクトルの大きさで角度を表すことができます。 よって、角速度・角加速度の表現にはもってこいなのです

慣性モーメント

慣性モーメントは回転において質量にあたる量です。

質量といえば重さと考えるかもしれませんが、理論的には力の影響の受けづらさを表す量です。 また、力が慣性を阻害する存在であることから、慣性への従いやすさを表す量であるともいえます。 少しラフな言い方をするならば「移動させにくさ」です。

慣性モーメントもまた、トルクの影響の受けづらさを表す量であり、回転の慣性への従いやすさを表す量です。 ラフな言い方をすると「回転させにくさ」になります。

慣性モーメントの値は回転軸周りの質量の分布によって決まります。 回転軸から遠い場所に高い質量があるほど回転させにくくなります。 しかし、慣性モーメントの値が回転軸の取り方によって変わってしまうとなると、 剛体の回転させにくさはどう表現すべきでしょうか。

そこで導入するのが「慣性モーメントテンソル(慣性テンソル)」と呼ばれる量です。 これは剛体のあらゆる回転平面に対する慣性モーメントをすべて表すことができます。 回転軸を決めて、慣性モーメントテンソルと演算することで、その回転軸周りの慣性モーメントを算出できるのです。

慣性モーメントテンソルは3x3の行列で表現することができ、 この行列に適切な回転行列を掛けると対角行列へと変換することができます。 この対角行列の対角成分のことを「主慣性モーメント」といい、回転後の座標系のことを「慣性主軸」といいます。

UnityではRigidbody.inertiaTensorが主慣性モーメントを、Rigidbody.inertiaTensorRotationが慣性主軸への回転を表します。 これらの値は質量とコライダの形状からそれらしい値がデフォルトで計算されますが、直接指定することもできます。

トルク

トルクは回転において力にあたる量です。 力のモーメントとも呼ばれます。

そもそも(n次)モーメントというのは距離のn乗に特定の物理量を掛けたものを指します。 トルクの場合は対象が力なので力のモーメントとなります。 慣性モーメントも英語ではmoment of inertiaなので直訳すれば慣性のモーメントです。

トルクは力によって生じる回転力、回転を加減速させる効果を表します。 トルクを与えると慣性モーメントの値に応じて角加速度が加えられます。

トルクは角速度や角加速度と同様にベクトルを使って表現します。

おわりに

移動方法の記事からコピペして書き換えた結果手抜き感がやばかったので、 おまけ程度に剛体の回転について簡単に書いてみました。 正直、自分もまともに勉強したわけではないので正しいかどうか怪しいのですが・・・。 もし何か間違っていたら連絡いただけると助かります。



執筆時のUnityのバージョン:2017.1.0p4

*1:ただし、AddTorqueの呼び出しではangularVelocityの値は即座に変化せず、内部で値が保持されたのちに物理演算時に反映されるようです

*2:密度や分布と言った方が正確かもしれません