2023年のコンパイル時レイトレーシング

C++

これは

qiita.com

の12/16の記事です。やばいもう16日が終わってしまう！

はじめに

太古の昔、あるC++プログラマ*1がtemplateの再帰と特殊化を使ってコンパイル時に計算ができることを「発見」*2*3したその日から、C++とコンパイル時計算は分かちがたく結びついています。

言語機能を本来意図されていなかった方法で活用する曲芸的な技巧だったコンパイル時計算は、その強力さを買われコミュニティを席巻し、C++11での constexpr の導入によりある意味で公式に認められたものとなりました。 導入当初は return 文一つだけしか持てなかった constexpr 関数は、それでも三項演算子による条件分岐や再帰によってコンパイル時計算を大いに盛り上げました。 そしてC++14で条件分岐、ループ、変数の書き換えが、C++17ではラムダが、C++20では仮想関数と動的メモリ確保が、とコンパイル時計算で利用できる言語機能は増え続けており、今では通常のプログラミング言語と遜色ないと言えるところまで来ています。

本記事では、まず constexpr の制限緩和の歴史を追いかけ、結果C++23でコンパイル時計算がどれほど簡単になったのかを、レイトレーシングを題材に見ていきます。 レイトレーシングはC++のコンパイル時計算では伝統ある（？）題材で、ボレロ村上さんの代表作の一つと言えるでしょう。 彼は芸術家にして日本で指折りのC++プログラマであり、またC++コミュニティ全体に絶大な影響を与えた constexpr 世界の開拓者でもありました。 レイトレーシングは、コンパイル時計算でどれほどのことができるのかを示す実例として、重い処理であることが想像でき、かつ視覚的インパクトを伴う形でその可能性を見せつける、いい題材だと思います。

今回はGitHubに実際に書いたものを上げてありますので、気になる方はそちらもご覧ください。 とはいえ、C++23が強力すぎて、ほとんど普通のコードになっていますが……。

constexprの発展史

最初はコンパイル時レイトレーシングのコードの解説だけで書こうと思ったのですが、書いてみたらあまりにも普通のコードになってしまったため、追加で constexpr の発展を時系列で見ることにしました。

ただ、見通しが甘すぎて長大になったので（甘すぎる）結構端折っています。 個別の要素に関してはもっとしっかりとした解説記事が別にあることが多いので、より詳細に知りたいという型は調べてみてください。

C++11

constexpr変数

constexprな変数は、コンパイル時に値が計算され確定し、実行時にそれを書き換えることはできなくなります。 これらの値はコンパイルの結果として即値としてエンコードされたり、.rodataセクションに置かれたりするでしょう。 constexpr指定できる変数はリテラル型に限られます。

constexpr関数

constexprな関数はコンパイル時に計算できる関数ですが、変数とは違って実行時に使うことも可能です。 結果をconstexpr変数にする場合はコンパイル時に実行され、このときは引数もコンパイル時定数である必要があります。

constexpr関数は変数を変更することはできませんでした。よって、計算結果は引数の参照経由などではなく、値として返す必要があります。そのため、voidを返すこともできませんでした。副作用がなく引数の変更もできないなら、voidを返す関数は何もできないからです。

さらに、constexpr関数の本体で許されるのは、型エイリアスやstatic_assertを除けば、return文だけでした。 よって、forやwhileなどの繰り返し計算やifによる分岐はできませんでした。 繰り返しは再帰呼び出しで、分岐は三項演算子cond ? then : elseで行う必要があります。

また、メンバ関数をconstexprにすると暗黙にconstになってしまいます。 C++11ではconstexpr関数内で変数を書き換えられなかったため、非constなメンバ関数は存在意義がなかったというところでしょうか。 この制限はC++14で、constexpr関数内での変数の更新が許可されると同時に緩和されました。

C++14

制限緩和

C++14ではconstexpr関数に大幅な制限緩和が行われました。具体的には、ifやswitchによる分岐、forやwhileによるループ、未初期化変数の宣言と変数の変更が許可され、ほとんど普通の関数になりました。

また、変数を変更できるようになったため、voidを返す関数も許可されました。引数を更新してvoidを返すような関数がconstexprになれるようになったからです。

さらに、変数を変更できるようになったので、constexprなメンバ関数が暗黙にconstになるという制限も取り払われました。constexpr文脈内で呼び出しても、クラスが変更される可能性が出てきたからです。

提案では、C++11の制限されたconstexpr関数で何かを書こうとするとかなりの面倒が生じること、D言語では当時既にこのような機能が存在していることなどが理由として挙げられています。*4 特に、標準ライブラリ関数をconstexpr対応させる際にも面倒が生じることが主張されています。例えば、イテレータのインクリメントは値の変更を（当然ながら）伴うので、constexprにできない、という点が指摘されています。

https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3597.html

<array>, <chrono>, <complex>, <initializer_list>, <tuple>, <utility>の一部

制限緩和と同時に、複数のライブラリ機能がconstexprに対応しています。

<array>のoperator[]とat()のうち、constなものがconstexprに対応しています。 これはC++11時点での提案だからでしょうが、C++14で変数の変更が可能になりconstexprメンバがconstにならなくなったというのに標準ライブラリはまだそれには対応していないというのは、少しやきもきするところではあります。

また、<chrono>の時刻演算がconstexprになっています。内部で使われている<ratio>はtemplate metaprogrammingによって実現されているため、コンパイル時計算との親和性はもとから高かったのかもしれません。

<complex>では、コンストラクタと比較演算子、それからreal(), imag()がconstexprになっています。 これも、C++11時点での提案だからでしょうが、real()やimag()のconstなバージョンだけがconstexprになっています。 演算などはまだサポートされていないようです。

<initializer_list>と<tuple>と<utility>、特に<utility>はstd::pairに対してconstexprが追加されています。 これらは、<complex>もそうですが、プレーンな構造体のように振る舞う型として対応されているように思えます。 それらのコンストラクタやメンバアクセスが、ユーザーが作るプレーンな構造体と同様にconstexprになってほしいというのは自然な要求でしょう。

std::tupleについてはarrayやcomplexとは異なり、std::getのmutable参照バージョンもconstexprになっています。 なんでそれを、std::arrayやstd::complexに分けてやれなかったんだ！　なんでなんでしょうね。提案されるタイミングの問題だったんでしょうか。

<utility>からはほかにも、std::moveとstd::forward（完全転送）がconstexprになっています。これも、constexpr文脈内での値の変更が可能になったためでしょう。

C++17

constexpr lambda

lambda式がconstexprで使えるようになりました。lambda式によって作られるクロージャオブジェクトをconstexpr文脈で呼べるようになり、またそれをconstexpr定数として扱えるようになりました。

lambdaにはキャプチャという機能がありますが、constexpr文脈で使うためにはキャプチャしたものは全てコンパイル時定数である必要があります……当たり前ですが。

とはいえC++17以前でも、constexprコンストラクタを持つクラスはconstexpr文脈で使えましたし、メンバ関数をconstexprにすることもできました。つまり、手動で対応するクラスを書けば、ラムダ式と同じことはできていたわけです。むしろ、このような不必要な回避策を取らなくて良いようににする、というのがmotivationであったようです。*5

if constexpr

少しconstexpr関数からは離れますが、関連する機能なので紹介しておきます。

これは、ifの後の条件がコンパイル時定数式である場合に使える機能で、選ばれなかった分岐を実体化しない（実行バイナリに含まない）というものです。

これによって、型特性に応じて使える最高効率のコードを生成するというような関数が簡単に書けるようになりました。具体例としては、IteratorのうちRandomAccess/ContiguousIteratorの場合は+Nを、そうでない場合はループを使ってインクリメントを続ける、というようなコードが簡単に書けます。 これ以前は、このような場合には別のタグ型を使ってオーバーロードによって切り替えていました。これは関数呼び出しの階層を増やしてしまうので、一回だけならまだしも、複数個の条件があると非常に大変になります。if constexprによってそれが一つの関数で、しかも普通のifの見た目で書けるのはかなりうれしいことでした。

ライブラリ対応

C++17で追加されたライブラリは、最初からconstexprに（できるものは）なっています。 例えばoptionalやvariant、string_viewなどです。もちろんいくらかの制限はあり、例えばoptionalは入っている値を変更することはできません。 anyですら、空のanyを作ること自体はconstexprになっていますね。流石に変更はできないようですが。

基本的にこの時点で、動的メモリ確保や入出力を含まないライブラリで、Cの遺産でないものはほとんどがconstexprになっているようです。ただし<algorithm>がconstexprになるのはC++20を待たなければなりません。

C++20

継承

これまで、仮想関数をconstexprにすることはできませんでした。 C++20では、constexprな仮想関数がconstexprな関数でオーバーライドされていた場合、定数式として使うことができます。 ただし、呼び出し側の基底クラスポインタとしてはどの関数が呼ばれる可能性もあるわけなので、オーバーライドしている関数は全てconstexprでなければなりません。

同時に、基底クラスポインタの先にどの派生クラスがあるのかあるのかを知るtype_idとdynamic_castも許可されました。

動的メモリ確保（new/delete、std::vectorとstd::string）

std::vectorとstd::stringがconstexpr対応し、またnewとdeleteもconstexprとなりました。

もちろん、格納する型はコンパイル時にコンストラクト・デストラクト可能である必要があります。 そしてnewのオーバーロードは許可されていません。

一番引っかかりやすい制限は、コンパイル時に確保された領域はコンパイル時に解放されなければならないということです。 灰は灰に。塵は塵に。コンパイル時に確保したstd::vectorを実行時に触れられるようにする方法はありません。 std::vectorやstd::stringを作ることはできますが、それはconstexpr関数の内部で使い切らなければならず、std::vectorを返す関数を作って、その戻り値をconstexprにすることはできません。

共用体

unionはメンバのうちの一つだけがアクティブな型です。再代入によってこのアクティブなメンバを切り替えるということが、constexprではこれ以前には許されていませんでした。

しかし、std::optionalやstd::stringはunionを使います。optionalはわかりやすいですが、std::stringはShort String Optimizationのために使っているはずです。std::stringをconstexpr対応させるときに、unionのアクティブメンバを切り替えられないことに気付いたようです。

これによって、std::optionalの代入関数やreset、std::variantへの代入などがconstexprになりました。

文法上の制限緩和

constexpr関数内では、たとえコンパイル時に呼んでいなかろうと、例えばインラインアセンブラやtry-catchを記述してはいけませんでした。constexpr内にtry {}が出てきた時点でコンパイルエラーです。 しかし、書くだけなら許してほしい（実行時に呼ばれることもあるのだから）という発想は自然です。

この提案で、try-catchとインラインアセンブラは、書くだけなら許されるということになりました。 これにより、if(std::is_constant_evaluated())を使うことで、実行時にはこれらの機能を使い、constexpr文脈で呼ばれた場合にはフォールバックを使う、というようなことが可能になりました。

ただし依然として、constexpr関数内で例外がthrowされるとコンパイルエラーとなります。将来的にはもしかするとcatchできるようになるかもしれませんが。

consteval

constexpr関数は、コンパイル時にも実行時にも呼ぶことのできる関数でした。

明らかにコンパイル時以外に呼ぶ必要のない巻数や、コンパイル時以外に呼ばれるべきではない関数も、constexprでは実行時にも呼べるようになります。それらの関数を実行時に呼んでしまうことを自動的に避けたい、あるいは、コメントよりもいい方法で読者に伝えたいという欲求は自然なものです。そしてconstevalはそれを実現します。

提案によると、どうやらstatic reflectionのことが念頭にあったようです。reflectionはコードの情報を使ったコードのことで、例えば今使っている変数のクラスが持っているメンバ変数の名前を全て得る関数であるとか、そういう機能のことです。 C++では静的なリフレクション、つまりコンパイル時にそれらの情報を得てコード生成を行うような機能が議論されており、そのような「コンパイル時以外に呼んではいけない関数」を明示する機能が欲しかった、というのがモチベーションの一つになっているようです。

is_constant_evaluated

std::is_constant_evaluated()は、コンパイル時に評価されている（つまり今がコンパイル中である）という場合はtrue、そうでない場合はfalseを返す関数です。 これを使ってifで分岐すれば、実行時にしか使えない機能、例えばインラインアセンブラによって狂気的にチューニングされたコードと、コンパイル時にも動く普通のコードを、一つのconstexpr関数内で同居させることができます。 実際、この機能のために、constexpr関数の文法上の制限緩和が進んだのでしょう。

ただし、if constexprとこれを組み合わせると、if constexprが評価されるのがコンパイル時であるせいで常にtrueになってしまうという大きな落とし穴があることに注意してください。 これはC++23で、if constevalによって解決されます。

<algorithm>と<numeric>、それに関係したライブラリ機能

基本的に渡されたイテレータ範囲で作業をする<algorithm>の多くの関数は、動的メモリ確保がなくとも（std::arrayがあるので）constexprで実装されていいように思えるものばかりでした。 ここで、それらの関数がconstexprになります。

とはいえ<algorithm>をconstexprにするために、いくらか追加の作業が行われていたようです。 たとえば、constexprで使えるIteratorの要件が定義されたり*6、swapがconstexpr対応したり*7とところどころが整備されていました。

また、<algorithm>の関数はものによってはコンパイラのイントリンシックや、もっとアグレッシブな場合はアセンブリによって高速化が図られていることがありました。これらの高速化をあきらめてconstexpr対応すると、当然ながら速度が落ちるのではないかということが懸念されていました。

提案では、速度低下が100%起きないと断言はできないものの、速度低下を回避しながらconstexprにできるだろうと述べています。その理由は、

  • 当時既にGCCの__builtin_memcmpはconstexprでも使えるので、コンパイラベンダーが頑張ることで他の組み込み関数、例えば__builtin_memmoveや__builtin_memsetも同様の処置ができるだろう。
  • アセンブリが使われている部分も、コンパイラのイントリンシックで置き換えることが可能だろう
  • 外部の関数が使われている部分も、それをインラインに展開するか、あるいはイントリンシックを使って再実装できるだろう
  • というか最近のコンパイラは賢いのでそういった最適化をヘビーにしている関数はそんなに多くないだろう

ということでした。

この提案によって、<algorithm>と<numeric>にある多くのアルゴリズムがconstexprになりました。 ここで見送られたのは、追加のメモリ領域を使う場合があるstable_sortなどのアルゴリズムと、まだconstexprになっていない乱数を使うアルゴリズムであるshuffleとsampleです。

https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2017/p0202r3.html

<complex>

かなり今更なのではないかと思うのですが、std::complexの演算がconstexpr対応しました。 isnanなどの<cmath>関数やコンパイラのビルトイン関数が使われていたことが、遅れた原因のようです。 というわけで、いくつかの簡単に実装できる数学関数しか使わない演算だけがconstexprになりました。四則演算と、normとconjです。

sqrtなどのより複雑な数学関数を使うstd::complex::absなどの関数は見送られました。

C++23

if consteval

ついこの間のアドベントカレンダー記事で紹介されていました。 一応説明をしておくと、これはconstexpr文脈にいるかどうかでコードを分けるための機能です。

if (std::is_constant_evaluated()) {/*1*/} else {/*2*/}を使ってコードを分岐させると、分岐は意図した通りにするのですが、実体にはどちらのパスも残ります。 もし/*1*/の中にconstevalなコードを持っていた場合、実際には通らないパスとはいえ、実行時に呼びだされる関数のコード内にconstevalが残るため、エラーになってしまいます。

ではこのような場合のために、通らないパスを実体から消し去るif constexprを、と思いたいところですが、if constexprはコンパイル時に処理されるため、is_constant_evaluatedを呼ぶと常にtrueになってしまいます。 この問題を避けるため、期待通りに動作するif constexpr (std::is_constant_evaluated())として、if constevalが導入されました。

ちなみにif constevalはC++20に向けた提案でしたが、締め切りを破っていたようです。それが原因なのかは知りませんが（多分そうでしょう）、これはC++23の機能になりました。

文法上の制限緩和

これまでも何度か（try-catchやインラインアセンブラ）ありましたが、C++23でもいくつかの文法機能がconstexpr関数内に存在することが許されるようになりました。

C++23では、constexpr関数内のラベルとgoto、リテラルでない型の変数、static, thread_localな変数の存在が許可されました。ただし評価してしまうとエラーになります。

https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2021/p2242r3.html

もう一つ、constexpr関数内でconstexprではない関数を呼んでいる場合、それが実際に評価されるまではエラーにならないようになりました。

これまで見てきたように、コア言語側でconstexpr実装が可能になったライブラリ機能がそれと同時に実際にconstexprになることは珍しいです。 大抵の場合は、1世代遅れてconstexpr指定されます。 もしconstexpr指定されていない関数を呼ぶconstexpr関数が、実際に呼ばれていないうちから警告を出していた場合、その警告を止めるために大量の機能テストマクロをconstexpr周辺に置くことになります。 これはまあ冗長でしょう。そういうわけで、constexprではないというエラーは、実際に評価されたときに出すようにするということになりました。

https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2022/p2448r2.html

簡単な数学関数

<cmath>の比較的単純な関数、つまりabsやround, ceil, floor、fmodとreminder、fminとfmax、それからfmaに、isfinite, isinf, isnanなどの判定関数、copysign, frexp, ldexp, scalbなどがconstexprになりました。

ただ、gcc以外のコンパイラはあまり対応していないようです（clang-18はダメでした）。

規格への変更がほぼconstexprの追加しかないので、このpaperを見ればどれが対応したのかはすぐにわかります。

https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2021/p0533r9.pdf

C++26

数学関数

<cmath>と<complex>に入っているほとんどの数学関数がconstexprになります。 これで、sqrtやsin, cosなどの普段使いの数学関数がconstexprになります。

https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2023/p1383r2.pdf

これで概ね、入出力を除くすべてのピースが埋まったように思えます。 シミュレーションやデータ解析などの計算で閉じているプログラムは、ほとんどがconstexprで簡単に動くようになったのではないでしょうか。

型消去

void*を通した型消去というテクニックがあります。 C++26では、これを行うためにconstexpr文脈内でのvoid*から別の型のポインタT*へのキャストが許可される予定です。 比較的知名度が低いテクニックのように思われるので、少し解説しておきます。

std::functionやstd::anyなどの「どんな型でも入る」クラスはこれを使って実装することができます。

std::anyのようにユーザーが取り出すべき型を知っていると思ってよい場合は、対象を格納する際にvoid*にキャストし、取り出す際にtemplate引数としてユーザーに型情報を教えてもらってそれにキャストして返します。std::any_cast<T>がそれに相当します。*8

std::functionのようにユーザーは取り出すべき型を知らない（忘れてよい）代わりにするべき操作が決まっている場合は、void*を適切にキャストして目的の操作をする関数への関数ポインタを覚えておくことでこれを実装できます。つまり、オブジェクト指向の説明でよくある Animal::speak() を例とすると、

class Cat
{
  std::string_view speak() const noexcept {return "meow";}
};

class Animal
{
  using func_type = std::string_view (*)(const void*);
  func_type speak_impl;
  const void* animal;

public:
  template<Speakable A>
  Animal(const A& a)
    : animal(std::addressof(a)),
      // 最初の型に戻してメンバを呼ぶ方法をここで作って覚える
      speak_impl([](const void* any) {
        return static_cast<const A*>(any)->speak();
      })
  {}

  std::string_view speak() const noexcept
  {
    return this->speak_impl(animal);
  }
};

のようにすれば、元の型に戻して対象のメンバ関数を呼ぶ関数を覚えておくことができ、ユーザーはもとの型Aを覚えておいて渡してやる必要がありません。

これだけでも便利ですが、他にも嬉しいことがあります。型消去をtemplateの代わりに使うことで、実体化されるコード量を削減できるという点です。templateの場合は型引数ごとに異なる実体が生成されてしまい、executableが大きくなりがちですが、型消去を使えば実体を一つになり、これを小さくすることができます。

これをconstexpr中でもできるようにしたいという提案がC++26に入る予定です。

https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2023/p2738r1.pdf

<algorithm>で前回スルーされていたもののいくつか

C++20時点で<algorithm>のほとんどの関数はconstexprだったのですが、一部の関数はそこに含まれていませんでした。 それら取り残された関数のうち、std::stable_sort, std::stable_partition, std::inplace_mergeがconstexprになる予定です。

これらはもともと、作業用のメモリを動的に確保したり、そこにplacement newするような操作を行うために後回しにされていました。*9 しかし既にconstexprで動的メモリ確保が可能になった以上、阻むものはなくなった形です。

https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2022/p2562r1.pdf

TMP

脱線ではありますが、template metaprogramming (TMP) についてもごく軽く説明しておきます。

TMPでは、再帰によって繰り返しを、（部分）特殊化のパターンマッチによって条件分岐を実現することで、templateの実体化を通して何らかの計算をすることを指します。

伝統に則って、素朴な素数判定を書いてみましょう。

namespace detail
{
template<std::uint32_t N, std::uint32_t I>
struct is_prime_impl
{
  constexpr static bool value = ((N % I) != 0) &&
    is_prime_impl<N, I-1>::value;
};

template<std::uint32_t N>
struct is_prime_impl<N, 1>
{
  constexpr static bool value = true;
};
} // detail

template<std::uint32_t N>
struct is_prime
{
  constexpr static bool value = detail::is_prime_impl<N, N-1>::value;
};
template<>
struct is_prime<1>
{
  constexpr static bool value = false;
};
template<>
struct is_prime<0>
{
  constexpr static bool value = false;
};

is_prime<N>::value には、Nが素数の場合trueが、そうでない場合falseが入ります。 どのような仕組みか見てみましょう。まず、このis_prime<N>::valueは以下のようになっています。

template<std::uint32_t N>
struct is_prime
{
  constexpr static bool value = detail::is_prime_impl<N, N-1>::value;
};

ここで、実際の値はdetail::is_prime_impl<N, N-1>::valueであることがわかります。そちらを見てみましょう。

template<std::uint32_t N, std::uint32_t I>
struct is_prime_impl
{
  constexpr static bool value = ((N % I) != 0) &&
    is_prime_impl<N, I-1>::value;
};

ここでは、NがIで割り切れるかどうかの値と、再帰的に展開されるis_prime_impl<N, I-1>::valueとの&&が取られています。 再帰部分では、Nはそのまま、Iは1ずつ減らして次の繰り返しに入っていることがわかります。 この再帰がI == 1のときに止まるべきなのはすぐにわかるでしょう。1は何でも割り切るので、常に(N % I) == 0になってしまいます。 再帰をI == 1で止めるのは、この部分です。

template<std::uint32_t N>
struct is_prime_impl<N, 1>
{
  constexpr static bool value = true;
};

これは「部分特殊化」と言って、第二template引数が1のときにマッチする型です。 これにより、I == 1となった場合はis_prime_implの実体はこちらになり、再帰展開は止まる場所を見つけます。 そしてここで変える値はvalue = trueなので、ここまでが全てtrue、つまりN % I != 0だった数は、true && true && .. && true = trueとなります。

また、これだとN == 1やN == 0だったときに困ります。なので、is_primeにはその場合のための特殊化を用意しておきます。

template<>
struct is_prime<1>
{
  constexpr static bool value = false;
};
template<>
struct is_prime<0>
{
  constexpr static bool value = false;
};

このように、再帰展開による繰り返しと、特殊化・部分特殊化によるパターンマッチを使うことで、テンプレートクラスの実体化を使って計算を行うことができます。

もちろん計算の対象は値だけでなく、using/typedefを使って型を出し分けることも可能です。 これと関数オーバーロードを組み合わせることで、特定の条件を満たす型（例：Iteratorのiterator_traits<Iterator>::iterator_categoryがRandomAccessIteratorであるもの）が来た場合にのみ関数の実装を分ける、というのは頻出のパターンになっています。 もちろん、現代ではこのようなことはConceptによって解決が可能となり、徐々に過去のものとなっていくでしょうが、過去の資産を紐解く際には便利な知識になるでしょう。

C++23でのコンパイル時レイトレーシング

一口にレイトレーシングと言っても、どこまでやるのかというのは難しい問題です。

今回は、非常に有名でありかつ既に多くの人が走ったことのあるレギュレーション、つまり"Ray tracing in one weekend"を走っていこうと思います。

https://raytracing.github.io/books/RayTracingInOneWeekend.html

この教材もC++を使ってはいますが、 何～故～か、実行時に走るコードしか書かれていないので、やる意味はあるでしょう。

とはいっても時間的な猶予と計算機リソースの問題があるので、最後のページの100個以上の球体が並ぶ画像まではいきません。中盤の、diffusive/metallic/dielectricなmaterialが適用された球体が並んでいるところまでとします。 実際、そこまでたどり着けば後は実行時間の問題でしかないと言えなくもないでしょう。*10

実際に出力された画像は以下の通りです。

コンパイル時計算で生成した画像。このクオリティだとコンパイル時間は大体30分程度。

画像ファイル出力

ファイルの入出力は流石にコンパイル時にはできません。 よって、ファイルのバイト列をコンパイル時に作っておいて、実行時にはそれがファイルに出力されるだけ、という状況にしておきたいところです。

std::vectorやstd::stringがconstexpr対応したとはいえ、コンパイル時に作ったstd::vectorを実行時に触ることはできません。 コンパイル時に作られたstd::vectorは、コンパイル完了と同時に解放されなければなりません。 ですので、ファイルはstd::array<std::uint8_t, N>にしておき、この先頭をconst char*として出力することにしましょう。

最初はppmのP3を使っており、数値を文字列に変換する関数を書いていましたが、どのコンパイラもconstexpr内で実行できる計算の回数に上限があるので、ここでそれを消費するのはまずいと考えP6を使うことにしました。 P6だとバイト単位でRGBが書き込まれるので、変換は必要ありません。 計算回数上限とメモリ使用量がほとんど唯一の縛りだったように思います。

std::arrayのサイズは関数内で決めて（もちろん定数ですが）、返り値はautoにしておきます。

constexpr auto
make_ppm(const std::array<pixel, IMAGE_SIZE_X * IMAGE_SIZE_Y>& image)
{
    constexpr auto header_string =
        "P6\n" STRINGIZE(IMAGE_SIZE_X) " " STRINGIZE(IMAGE_SIZE_Y) "\n255\n";
    constexpr auto header_size = std::string_view(header_string).size();
    constexpr auto pixel_size = 3; // [byte]
    constexpr auto ppm_size = header_size +
        (IMAGE_SIZE_X * IMAGE_SIZE_Y * pixel_size);

    std::array<std::uint8_t, ppm_size> retval;
    for(std::size_t i=0; i<header_size; ++i)
    {
        retval[i] = static_cast<std::uint8_t>(header_string[i]);
    }
    std::size_t idx = header_size;
    for(std::size_t pxl=0; pxl<IMAGE_SIZE_Y * IMAGE_SIZE_X; ++pxl)
    {
        const auto [r, g, b] = image[pxl];
        retval[idx++] = r;
        retval[idx++] = g;
        retval[idx++] = b;
    }
    return retval;
}

これが成功しているかどうかは、出てきたバイナリを見るとわかります。 xxdで見てヘッダのP6を検索してみると……。

実行バイナリをxxdで見ているところ。.rodataセクションにPPM画像のヘッダが入っており、おそらく画像本体と思しきバイト列が続くことが見て取れる。

となり、バイナリに実際に画像データが埋め込まれていることがわかります。

数学関数

レイトレーシングではベクトルの距離を計算する必要があるため、sqrtは必須です。 1 / \sqrt{x}を計算するrsqrtがあってもよいでしょう。

sqrtのような数学関数がconstexprになるのはC++26からなので、C++23では使えません。これに関しては自力でどうにかする必要があります。 今回はC++23でも使える<cmath>関数であるldexpとfrexpを使って浮動小数点数を指数部と仮数部に分解し、仮数部を256分割してそれぞれの区間を2次関数で補完することにしました。

別のコードでパラメータを生成しておけば、2次関数の計算は本当にやるだけです。

constexpr inline std::size_t sqrt_N = 256;
constexpr inline double sqrt_dx = (1.0 - 0.5) / sqrt_N;
constexpr inline double sqrt_a[sqrt_N] = {
#include "parameters/sqrt_a.dat"
};
constexpr inline double sqrt_b[sqrt_N] = {
#include "parameters/sqrt_b.dat"
};
constexpr inline double sqrt_c[sqrt_N] = {
#include "parameters/sqrt_c.dat"
};

constexpr double sqrt_interp(const double x)
{
    const std::size_t i = (x - 0.5) / sqrt_dx;
    const auto a = sqrt_a[i];
    const auto b = sqrt_b[i];
    const auto c = sqrt_c[i];
    return (a * x + b) * x + c;
}

sqrtの場合、指数部は単に半分にするだけで計算できます。奇数だった場合は、最後に余った1に相当する \sqrt{2}をかける必要があることに注意しましょう。また、指数が負だった場合には逆に割らなければなりません。 infやnanの場合についても、C++23時点で既にisinfやisnanがconstexprなので特に問題はありませんでした。

現実的には、clangがfrexpとldexpに対応していないので、これを実装する必要があります。 とはいえstd::bit_castがconstexprなので、doubleをstd::uint64_tに変換したのちビット演算で指数部と仮数部を取り出して、仮数部を[1, 2)から[0.5, 1)に変換すればよいだけなので、自作は容易です。

指数部が負になっている場合などの扱いをかなり愚直にやったので長くなってしまっていますが、以下のようになります。

constexpr double sqrt(const double x)
{
    if(x == 0.0 || x == -0.0)
    {
        return 0.0;
    }
    if(std::isnan(x))
    {
        return std::numeric_limits<double>::quiet_NaN();
    }
    if(std::isinf(x))
    {
        return std::numeric_limits<double>::infinity();
    }

    int ex=0;
    const double norm = std::frexp(x, &ex);

    const int zex = ex / 2;
    if(ex >= 0)
    {
        if(ex % 2 == 1)
        {
            const double z = sqrt_interp(norm) * std::numbers::sqrt2;
            return std::ldexp(z, zex);
        }
        else
        {
            const double z = sqrt_interp(norm);
            return std::ldexp(z, zex);
        }
    }
    else
    {
        if(ex % 2 == -1)
        {
            const double z = sqrt_interp(norm) / std::numbers::sqrt2;
            return std::ldexp(z, zex);
        }
        else
        {
            const double z = sqrt_interp(norm);
            return std::ldexp(z, zex);
        }
    }
}

乱数

<random>はconstexpr対応していないので、乱数に関しては自分で作る必要があります。 そこまで高品質のものを求めているわけではないので、簡単にxorshiftでよいでしょう。

struct xorshift64
{
    std::uint64_t state;
};

constexpr xorshift64 next(const xorshift64& curr) noexcept
{
    std::uint64_t x = curr.state;
    x ^= x << 13;
    x ^= x >>  7;
    x ^= x << 17;
    return xorshift64{x};
}

どちらかというと問題はここよりも、これをdouble上の乱数にすることです。 とりあえず[0, 1)の一様乱数が作れるとほとんどのことが何とかなるので、これを作ります。

今回は簡単に、指数部を0に固定して仮数部を乱数にすることで[1,2)の乱数を作り、ここから1を引いて[0, 1)にすることにします。 これだと得られる値が仮数部の最小単位―doubleだと2^-52―刻みになるので、本来指数部を小さくすると表現できるはずのより小さな値が表現できないという問題がありますが、今回は特にそこは必要がないので、ここで切り上げます。

std::bit_castは最初からconstexprなので、uint64_tで作ったビット列をそのままdoubleに変換することができます。

constexpr std::pair<double, xorshift64> uniform_12(const xorshift64 rng) noexcept
{
    static_assert(std::numeric_limits<double>::is_iec559, "assuming IEEE 754");

    const auto nrng = next(rng);

    const std::uint64_t mantissa_mask   = (1ull << 52) - 1;
    const std::uint64_t random_mantissa = nrng.state & mantissa_mask;
    const std::uint64_t exponent_bias   = 1023ull << 52;
    const std::uint64_t float64         = exponent_bias + random_mantissa;

    return std::make_pair(std::bit_cast<double>(float64), nrng);
}
constexpr std::pair<double, xorshift64> uniform_01(const xorshift64 rng) noexcept
{
    const auto [d, r] = uniform_12(rng);
    return std::make_pair(d - 1.0, r);
}

勢い余ってxorshift64自体を毎回更新していますが、別に変更しながらでもよかったですね。まあ特に変わりはありませんが……。

浮動小数点数については以前に記事を書いているので、どういう表現になっているのか知らないという人はこちらもご覧ください。

in-neuro.hatenablog.com

in-neuro.hatenablog.com

オブジェクトのマテリアル

レイトレーシングにおいてオブジェクトには、位置と形状の他に、それに光が反射したときに反射光の色がどう変化するかの値と、どの方向に光が反射するかの情報が入っています。 色は単純にRGBがそれぞれどの程度減るかという情報で管理されていますが、光の反射は分布も計算方法もかなり異なるので、異なる関数を呼び分ける必要が生じます。

オブジェクトによって同じ名前の関数の振る舞いを変える方法は大きく分けて二つあります。std::variantと継承です。 継承もC++23ではconstexprで使えるのですが、今回は単純な好みでstd::variantを使いました。

std::visitを使ってconstexpr lambdaでオーバーロードされているscatter関数を呼べば、マテリアルに応じたscatter関数が呼ばれます。

struct object
{
    color attenuation;
    std::variant<diffusive, metallic, dielectric> material;
    std::variant<sphere> shape;
};

template<std::size_t N>
constexpr std::pair<color, xorshift64>
ray_color(const ray& r, const world<N>& w, const xorshift64 rng, const std::size_t depth)
{
// ...
        const auto [nray, nrng] = std::visit([&r, &collision, rng](const auto& mat) {
                return scatter(r, collision, mat, rng);
            }, obj.material);
// ...
}

実行

g++-13とclang++-18を使って実行しています。 コンパイル時計算をするにあたっては、どちらも一長一短あります。

例えば、gccとclangはどちらもコンパイル時間が長くなりすぎないようにコンパイル時計算で許す計算の回数を絞っているのですが、この回数の上限が異なります。gccは-fconstexpr-ops-limit=Nで、clangは-fconstexpr-steps=Nでこの上限を変更することができます。 できるのですが、gccは2^32-1を超えても問題ないのに対して、clangは232以上の値を指定することができません。 よって、clangの公式バイナリでのコンパイル時計算はgccと比べると著しく制限されています。 ただし、この制限解除には多くの人が取り組んでおり、以下のような変更を加えることでこの上限を解除したclangをビルドすることができます。*11

qiita.com komorinfo.com

また、clangではldexpやfrexpがconstexprになっていないなど、マイナーな部分とはいえgccと比べて少し対応が進んでいない部分があります。

ではgccを使っておけば何の問題もないのかというと、そうでもありません。 確かにgccではconstexprでの評価回数制限がかなり緩いのですが、代わりにすさまじいまでのメモリを消費します。 内部的に定数式の値を全て覚えているのでしょうか？　秒間数百MBのペースでメモリ使用量が増えていき、今回の最後の画像を生成する時は128GBを突破してしまいました。 これ以上のメモリを積んだマシンは用意できなかったので、gccを使っても無制限にコンパイル時計算ができるということはありません。

というわけで、おそらく最善の方法はパッチを当てて-fconstexpr-stepsの上限を解除したclangを使うことだと結論付けられます。

実際のコード

コードは以下で公開しています。

github.com

コンパイルに時間がかかりすぎるので、CIを設定する気になれませんでした……。

おわりに

ちょっとコードがあまりにも普通になったのでconstexprで縦軸でまとめるか、と思ったら非常に大変なことになってしまい（見通しが甘すぎる）、クソ長くなった上に死ぬほど遅れました。 とはいえ、たまに自分でもどこまでがconstexprで使えるかを覚えていないことがよくあったので、いい機会だったと思います。 実際には（特にライブラリ機能に）見落としがあると思いますし、ひょっとすると誤解してるものもあるかもしれません。その場合はあとで更新するかもしれません（した場合は一番上に書いておこうと思います）。

書いてみると、C++23ではコンパイル時計算はほとんどもう普通のコードになっていることがよくわかります。 事実、入出力と乱数や数学関数などがまだ対応していない、ということを念頭に置いて普通のコードを書くと、それがそのまま動くというような状態でした。

一番最初に一回計算しておけば済むことを先に計算する、というような目的で使っている場合、数学関数を除けばほとんど障害になるものはないでしょう。 制限だらけだったころのことを思うと少し……なんというか……面白くなくなったという気はしますが、実際に書かないといけないコードを書いているときは楽に書けた方が良いに決まっているので。 便利な時代になったと思います。

*1:Erwin Unruh, http://www.erwin-unruh.de/

*2:Item 48: Be aware of template metaprogramming. Effective C++: 55 Specific Ways to Improve Your Programs and Designs, Scott Meyers

*3:History of TMP, C++ Programming/Templates/Template Meta-Programming, https://en.wikibooks.org/wiki/C%2B%2B_Programming/Templates/Template_Meta-Programming#History_of_TMP

*4:https://dlang.org/spec/function.html#interpretation

*5:https://isocpp.org/files/papers/N4487.pdf

*6:https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2018/p0858r0.html

*7:https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2018/p0879r0.html

*8:格納する際に簡単なのは動的メモリ確保をしたのちそこへのポインタをvoid*にすることでしょう。実際の std::any は、Small String Optimizationに似た最適化によって動的メモリ確保を回避しています

*9: III C. Analysis of existing algorithm implementations. https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2018/p0879r0.html

*10:とはいえそこはコンパイル時計算では結構本質的に重要な問題ではあるのですが……。

*11:今回やろうと思っていたのですが時間がとれずできませんでした……。