Syoyo Fujita's Blog

raytracing monte carlo

Month: October, 2007

Youngest Prof. nyaxt recieves “incredible creator” title

[En]

Prof. nyaxt is honored as “incredible creator” by Exploratory Software Project sponsored by IPA [1],
with his physically-based and unified renerer “nytr”,

http://nyaxtstep.com/projects/nytr

In Exploratory Software Project, very few people is titled as valuable “incredible creator” for their extraordinary talent.
And more surprisingly, nyaxt is the youngest “incredible creater” ever since Exploratory Software Project started in 2000.
Yah!, he is definitely genius and newtype(if you don’t know about the word `newtype’, wikipeding it [2])
in the world of renderer writers.

[Ja]

nyaxt 先生が、氏の物理ベースなユニファイドレンダラでスーパークリエイターに認定されています。
が、さらに驚いたのはこれ。

http://it.nikkei.co.jp/business/news/index.aspx?n=MMIT0z000029102007

情報処理推進機構(IPA、東京・文京)が優秀なプログラマーを支援する「未踏ソフトウェア創造事業」の「天才プログラマー/スーパークリエータ」に、千葉大学理学部2年の nyaxt さん(18)を認定したというのがあった。2000年度に始まった同制度で最年少の認定だという。

し、しらなかった、 nyaxt 先生が最年少のタイトルホルダーだったなんて…
いやー、これはきっとニュータイプに違いない。

最近ではわずか 2 ヶ月でパストレを書いちゃうひともいるし… [3]
最近のレンダラ野郎はすごいやつばっかりです。
これからのレンダラ野郎はもしかしてみんなニュータイプなのでは?と思ってしまいます。

少し残念なのは、この記事を書いた人や PM が、(大域照明)レンダラのことが分かっていない
(正当に評価できていない)こと。
この程度の評価は低すぎで、もっと正当に高く評価されてしかるべきだと思います。

「本当にコワイ奴は 下から来るんだ」

私は今、まさにリアルに○○○の碁の○田六段の心境です。

オレもうかうかしてらんねー、レンダラの巨塔の頂点を目指していればいいかと思いましたが、
これからは下から這い上がってくるヤツにも本気で注意しないとです。
オレはニュータイプじゃないから、ニュータイプに来られたらマジ終わる。
まぁすでに何人かには抜かれちゃいましたけど…

そうでなくても、世界に目を向ければ 2050 年には世界一の IT & 経済大国になるインドの
インド人を初め、 EU 人や、南米などの新興国人も脅威の存在です。
サブプライムショック -> メリケン景気減衰により、新興国の発展が予定よりも早く加速してきています。
これからものすごい勢いで台頭してくることでしょう(EU はドル安ユーロ高による相対的要因)。

[1] IPA is (possibly) biggest independent administrative corporation for IT technology in Japan.
http://www.ipa.go.jp/index-e.html

[2] http://en.wikipedia.org/wiki/Newtype

[3] ktkray
http://d.hatena.ne.jp/ktkray/

Advertisements

Dive into the tumblr world…

… tumbr が面白いのでしばらく tumbr の世界で過ごすことにします.
ついでに、tumbr サービスは 11/1 になにかあるみたい。

http://blog.davidville.com/

まあ tumbr ばっかやっていても情報創造には寄与しないと思うので、
いずれ戻ってこないといけませんが。

I’m now Tumblelogging

[En]

I just started my Tumblr at

http://lucillerender.tumblr.com/

Scraps, small infos, quotes and more on small thing about developing my lucille renderer is now found there, not at this blog.

[Ja]

ただでさえ内容の品質が低いこの開発日記ですが、
最近は数行でリンクを紹介するだけの情報創造 [1] にほとんど貢献しない
しょーもないエントリも増えてしまい、もっと品質を落としてしまっています。
読んでいただいている方に申し訳なく思います。

とはいえ、自分の備忘録も兼ねて、いくらかは記録をとっておきたい。
日記の品質を落とさず必要そうなネタはとっておきたい。
どうしたらよいかと思っていたら、ちょうど Tumblr(たんぶらー) という良さげなツールを見つけました。
使い方は [2, 3] などを参照。

というわけで、これからは小さなネタ、ニュースとか web スクラップとかは Tumblr のほうに tumblelogging
して、それをかもした(?)ものをこちらの日記に書いていきたいと思います。

[1] 貨幣の信用創造の情報版
[2] これから15分でtumblrを始めるための資料
http://blog.overlasting.net/2007-05-12-2.html
[3] Tumblrの魅力:かつてないほど「Webをやる」ことに惑溺できるツール
http://nakanohajime.wordpress.com/2007/10/27/tumblr/

統計科学のすすめ[その2]

susemi07-11.jpg

統計科学のすすめ[その2]
http://www.nippyo.co.jp/maga_susemi/ss0711.htm

[En]

“Math seminer” is a math magazine published in Japan.
In the issue of this month, the manazine features cutting edge statistical science,
which is much valuable for thinking MCMC based rendering for the next of MLT.

[Ja]

たまたま本屋で見かけて、統計科学について特集していたのでおぉっと思いつつ手に取ってみると、
予想通り? 伊庭先生がトップで解説していたので早速購入。
(「その1」 は残念ながら知らなかった。バックナンバーを図書館で探してみよう)

「これからは統計学も個性を重視しなければいけない」

うーむ、なるほど。基本的にはグローバルなパラメータでモデリングするんじゃなくて、
ローカルなパラメータでモデリングできるようにしよう、でもあまりローカルすぎるのも
だめだよ、と私は理解してみました。

spherical harmonics は global support だけど、wavelet は local support だよ、
でもじゃあ wavelet でいいんじゃないの? かというと、wavelet には周波数空間での効率的な
rotation が出来ないからやっぱ spherical harmonics がいいときもあるのです、
みたいな感じでしょうか。

たんなる統計学(統計的平均など)がまったく世の中の現象で役に立たないことは、
特に金融市場において、たとえば LTCM 危機やサブプライム危機で証明されています。
また、plain vanilla な black and scholes 方程式なども、結局は理想式なわけで、
現実では volatility smile などにより市場価格をうまくトラッキングできないわけです。

レンダラの世界も、基本は非積分関数が discontinuity で unbound なので、
ふつうの統計的な手法はやっぱりあまり期待できないと感じています。

じゃあどんなのがいいかというと、具体的にはまだあまり思い浮かびませんが、
アイデアとしては、LTCM 危機などのように、
「統計的にあり得ないことが実はとても起こりやすい」ことを
前提としたモデリング手法みたいなのがあるといいのかなぁと。
まだよく知りませんが、金融系では jump diffusion とかというモデリング方法がそれにあたるのかもしれない。

あとは最近は Adaptive MCMC [1] とかが出てきているので、
これも役に立ちそうでしょうか。

おまけ

こんどはメリルが四半期で 1 兆円近い損失ですか [2] …

… そんななか、日本の農林中央金庫がやってくれました![3]
さすが農中!俺たちにできないことを平然とやってのける!
そこにしびれる憧れるゥ!

これがうまくいくと、がぜん面白くなってきますね。
いやー、やっぱクオンツの道を歩むべきだったのかなぁ…

[1] MCMC Preprint service.
http://www.statslab.cam.ac.uk/~mcmc/pages/list.html

[2] 米メリル:7―9月期は2001年以来の赤字、79億ドルの評価損計上
http://www.bloomberg.com/apps/news?pid=90003009&sid=aiaTbgVpev.M

[3] 農中:欧米のABSに3兆円追加投資へ-サブプライム余波で妙味(3)
http://www.bloomberg.com/apps/news?pid=90003017&sid=aRex_iHUaZ1g

独りで始める Concrete and Specific Programming(CSP)

独りで始める Concrete and Specific Programming(CSP)
http://madscientist.jp/~ikegami/diary/20070822.html#p01

このエントリでは、開発側の夢を実現するためのプログラミング技術 Concrete and Specific Programming を提唱します。
Concrete and Specific Programming では、最初に「夢」を書くことを提唱します。次に、設計書を書き、次に、仕様に基づいたテストスイーツを書き、最後に実装します。

Haskell の QuickCheck 関連でヒットしたページ。
私は今現状の lucille を gappa や coq, SPIN などの Formal verification や model verification を使って
一から書き直したい要求にかられているのですが、
この CSP の考え方がとても共感できます。
CSP のやり方と formal verification の考えが似ている気がするので、
少しテストについて考えてみたいと思います。


テストについて

いままではコードに対するテストでいいかなぁとか思っていたのですが、
gappa を知ってから考えが変わりました。
レンダラにおいては、アルゴリズムの検証が非常に重要になると。

そして、コードは仕様やアルゴリズムの実態であると考え、
SCP のように仕様に対するテストを書く、アルゴリズムに対するテスト(証明)を書く、
というものが必要である、と。

たとえば、

半球上に点を一様分布生成する

コードを実装したいと思います。

いままでの私なら、

– 本とか参考にしてアルゴリズムを実装して
– ほんとうにそのコードの結果が、一様分布する点群を生成しているのかテストコードを書く

という手順を取っていました(formal verification を知らなかったから)。
一様分布しているかどうかのテストはたとえば、
コードを通って出てきた結果を gnuplot で目視したり、
二点の距離を測るプログラムを書いてとあるしきい値以下であることを確かめる、
になります。

ただ、このようなやりかたってどうしても抜けがでる。
テストがすべての条件をカバーできているのかよくわからないし、
そもそもアルゴリズムが間違っていて、たまたまテストでの条件ではうまく振る舞っているだけなのかも
してません。

それにそもそもこういう感じでコードに対するテストを書くというのがめんどくさい。
めんどくさい上に、下手すると brute force にテストする場合もあるかもしれません。
そんなめんどくさいのに時間を取られていたらインド人に追い抜かれちゃうよ。
日本人はこれからどんどん知的生産性を上げていかないと、マジでこの世の中生き残れません。

アルゴリズムをうんうん考えて、間違えてないかチェックして、チェック OK ならあとは
コードに落として完成、とやってアルゴリズムの思考に集中できないものか。

そんなわけで、SCP や formal verification 系では、アルゴリズムに対してテスト(証明)を書く。

– アルゴリズムを記述して
– そのアルゴリズムが正しいかどうかテスト(証明)する
– テストが OK(証明ができれば)したら、コードに落とす.

コードはテスト(証明)されているので、基本的にバグなしでありきちんとアルゴリズムを実現できているものになる。

つまりテストをもっと上のレベルで行うということですね。

SCP ではテスト手法としては、QuickCheck などのランダムテストを上げています。

いろいろ調べてみると formal verification や model checking などのキーワードで、
形式的検証という分野があることがわかりました。
除算のバグで 5 億ドルの損失を出した会社は、以降、プロセッサの設計には
この形式的検証を重視するようになったそうです[2]。

ソフトウェアプログラミングの方面でこれら形式的検証の
具体的なツールとしては、gappa(浮動小数点計算検証)、SPIN(モデル検証)、
Coq(証明支援)、Isabelle(自動証明)、caduceus(C のコードを検証) などがあるみたい。
これらを仕様やアルゴリズムを検証するためのツールを使うとより効率的な気がする。

SPIN は NASA のシステムでも使われていたりして、また歴史も古く実績はあるようだ [3]。
通信プロトコルの検証にも使えるみたいなので、並列レンダリングの処理に使えそう。

とはいえ、レンダラ作成においてアルゴリズムが重要であるとはいえ、
仕様•アルゴリズムレベルでテストやって OK が出たら、あとはコードに落とすだけ
という仕組みで全部物ごとが運ぶかというと、そうとも言い切れません。

コードレベルでの最適化や省メモリ処理など実装中心の問題もでてきます。
レンダラ作成においてはアルゴリズムもコードの速度も重視すべき点であるので、
ある程度どこかでトレードオフをすることになるでしょう。

すくなくとも C での SIMD 最適化コーディングは MUDA で解決できるようにして、
SIMD intrinsics を使って C コードを最適化する、
という非生産的なコーディングステップはなくしたいですね。
(SIMD プログラミングをしっかり学びたい、というフェーズの場合は除く)

[1] Gappa : automatic proof generation of arithmetic properties
http://lucille.atso-net.jp/blog/?p=338

[2] 自動定理証明(wikipedia)
http://ja.wikipedia.org/wiki/自動定理証明

[3] FM2003 報告
http://www.shinsahara.com/www/sigfm/material/FM2003%95%F1%8D%90.pdf
「NASAってSPIN実際に使ってるんか?」
「全宇宙システムで使っている」

floating point computation, part 0: 数値計算と誤差の話

渡部 善隆:
数値計算と誤差の話〜浮動小数点演算はどれくらい信用できるか〜,
九州大学大型計算機センター広報, Vol.27, No.5 (1994), pp.556-580.
http://yebisu.cc.kyushu-u.ac.jp/~watanabe/RESERCH/MANUSCRIPT/KOHO/FLOATING/intro.html

最近浮動小数点計算の資料を読みあさっていて、一番印象に残った文がこれ。

数値計算の品質保証に関する研究は、
ヨーロッパ、特にドイツで精力的に行なわれて います。
もし、読者のどなたかがヨーロッパの学会で、
例えばある非線形方程式の数値計算などの研究発表をされる場合、
精度保証付のソフトウェアを使って誤差をきちんと評価した結果を出してないと、
ましてや単精度の計算などだったら ブイブイ突っ込みが来る可能性のある分科会もある
ので、
特にドイツに行かれる方はご注意下さい。

なんと胸を突き刺すようなお言葉…
この文で言われているように、浮動小数点計算の文献を多く当たって感じるのは、
単純に計算が早ければよいわけではないということ。
そして、単精度ではかなり注意しなければ精度に不安がありそうだということ。

グラフィックス、とくにレンダリング、の研究発表において、
数値計算の誤差について気に留めるやつなんて見たことが無い気がする。
(ジオメトリやシミュレーション系ではあるかと思いますが)

たとえば GPGPU で、
「GPU で CPU の浮動小数点計算が n 倍早くなりました!」というたぐいのもので、
ただし GPU の実装では誤差が x ulp になります、それに対して CPU では y ulp ですみたいな解析はやっているんでしょうかね。
(実行してみたら誤差が 1.0e-4 ですとかそんなのは当たり前として)
CPU は倍精度で計算していて、GPU は単精度だった、なんてことも無きにしもあらずな世界ですからね。

とはいえ、最近ではしっかりしたヤツも現れているようで、
[3] 見たいに GPU での量子計算で和を取る項の計算が長くなるので精度を上げるために
Kahan summation を使いました、というものもあります。
この論文では Kahan summation を使っていますが、
GPU の馬力を生かして Ogita, Rump, Oishi の最強の誤差なし総和計算アルゴリズム [4] を
使ってもいい気がします。

オフラインレンダラにおいては、浮動小数点計算がいい加減で、
ロケットが爆発したり 、5 億ドルの損失を出して 4 半期の利益の半分を吹っ飛ばした [1]、
という事態になることはないですが、
何時までたってもレンダリング画面の輝点や黒点の原因が特定できずに
デバッグに莫大な時間を費やしてしまうということはあり得ます。
(lucille の今の実装がまさにそれ. そして分かったからといっても、
ではどう対処するとどれくらい誤差を抑えられるのかの理論がよくわかっていないので、
対処のためのモデル式が立てられない)

数年前、とあるインハウスプロダクションレンダラを書いているレンダラ神に聞いた話では、
かなり浮動小数点の計算には注意を払っているとのことでした。
プロダクション用ですから、輝点や黒点とかが出たら致命的ですからね。

当時はでも、ふーんそうなんだ、と思うくらいで
あまり気に止めていませんでした。
しかし別のところで、 CPU が違うと
浮動小数点計算の違いでレンダリング結果が異なるんですよ、
とか聞いたりして、浮動小数点計算の問題は根が深くきちんと対応しないといけないな、
と思うようになりました。

浮動小数点計算の問題をないがしろにしてきたツケが今になって自分に大きく響いてきています。

もっと浮動小数点計算のことを知っておく必要があります。
というわけでレンダラやろうのための浮動小数点計算をシリーズもので企画中。
いくつかのメモは google note に書きとめて、[2] に残してあります。

おまけ: 我がドイツ産業銀行のォ〜損失額はァ〜〜、世界一ィイイイイイ!!!!

ヨーロッパでは数値計算の品質保証に厳しいということなのに、
ここ最近のヨーロッパ系銀行のサブプライムローンによる損失はなんなんでしょうね。
IKB の SIV がデフォルトなんて起こしているし… これってかなりやばいでしょ[5,6]
もしかして 2 兆円以上の損失!?
ドイツ銀行も 4 半期で 4, 000 億近い損失だしてるし…
もう桁がちがうよね、LTCM の時代は遠くになりけり…

さすがおいもの国!俺たちにやってのけられないことを平然とやってのける!
そこにしびれる憧れ…るわけないか。

リスク管理と数値計算の品質保証の分野の繋がりはあまりないということなんでしょうか。

[1] http://wiredvision.jp/archives/200511/2005111605.html
http://findarticles.com/p/articles/mi_m0EKF/is_n2049_v41/ai_16811996

[2] 浮動小数点計算について
http://www.google.com/notebook/public/13822895759654633739/BDQbSIgoQ37K-u9Ai

[3] Quantum Monte Carlo on Graphical Processing Units
Amos Anderson, William A. Goddard III, Peter Schröder,
Computer Physics Communications, 177(3), pp 298-306
http://www.multires.caltech.edu/pubs/

[4] Accurate Sum and Dot Product
Takeshi Ogita, Siegfried M. Rump, Shin’ichi Oishi
http://citeseer.ist.psu.edu/711647.html

[5] IKB傘下のSIVラインブリッジ、条項発動で全債務の即時返還必要
http://www.bloomberg.com/apps/news?pid=90003017&sid=atYLIv.muKlw&refer=jp_japan

[6] 英チェーンCと独IKB銀運営のSIVがデフォルト-資産価値下落で
http://www.bloomberg.com/apps/news?pid=90003009&sid=alMC9xGW7Brg&refer=jp_home

luxrender, a GPL unbiased renderer

http://www.luxrender2.org/

LuxRender is a new, open-source, free software rendering system for physically correct, unbiased image synthesis.
Rendering with LuxRender means simulating light according to physical equations, this produces realistic photographic quality images.
It’s an authorized fork of the PBRT project, focusing on production rendering and artistic efficiency.

PBRT based unbiased renderer. Cross platform, multithreading, MLT, SSE2 optimization and etc… very nice!
(but one disappointing thing. it is coded in obsolete C++ with boost combination.)

Work in progress: Initial Haskell version of MUDA

[En]

Finally, I’ve finished initial Haskell reimplementation of MUDA language[1,2, 3] compiler.

http://lucille.svn.sourceforge.net/viewvc/lucille/angelina/haskellmuda/

Here is an example.


// input.mu
// Input MUDA code
vec add_func(vec a, vec b)
{
        return a + b;
}

From this input, MUDA compiler(with x86 SSE CodeGen) generates the following SIMD code.


$ mudah input.mu

//
// The following code is generated by MUDA compiler
//
#include 
#ifdef __64bit__
    #error "64bit env is not yet supported"
#else
#endif

#ifdef __GNUC__
static inline void *muda_aligned_addr16(void *addr) {
    return (void *)((((unsigned int)addr) + 15UL) & ~(15UL));
}
#endif

#ifdef __GNUC__
    #define MUDA_ATTRIB_ALIGN __attribute__((aligned(16)))
    #define MUDA_DECL_ALIGN
#elif defined(_MSC_VER)
    #define MUDA_ATTRIB_ALIGN
    #define MUDA_DECL_ALIGN __declspec(align(16))
#else
    #error "Sorry, MUDA doesn't support your compiler
#endif

static inline __m128 muda_sel_ps( __m128 a, __m128 b, __m128 mask )
{
    b = _mm_and_ps( b, mask );
    a = _mm_andnot_ps( mask, a);
    return _mm_or_ps( a, b );
}

__m128 add_func (__m128 a, __m128 b)
{
    const __m128 t_vec1  =  _mm_add_ps( a ,  b  ) ;
    return t_vec1 ;
}

And more, another example. SIMD ray-triangle intersection code written in MUDA.


// isect.mu
//
// SIMD cross product
static inline vec cross4(vec a, vec b, vec c, vec d)
{
        return ((a * b) - (c * d));
}

// SIMD dot product
static inline vec dot4(vec ax, vec ay, vec az, vec bx, vec by, vec bz)
{
        return (ax * bx + ay * by + az * bz);
}

//
// 1 ray - 4 triangles intersection.
// Assume that, for ray data, same value is copied over 4 elements.
//
vec
isect(vec rox, vec roy, vec roz,
      vec rdx, vec rdy, vec rdz,
      vec v0x, vec v0y, vec v0z,
      vec e1x, vec e1y, vec e1z,
      vec e2x, vec e2y, vec e2z,
      out vec outT, out vec outU, out vec outV, vec prevT)
{
        vec px, py, pz;

        // p = d x e2
        px = cross4(e2z, e2y, rdy, rdz);
        py = cross4(e2x, e2z, rdz, rdx);
        pz = cross4(e2y, e2x, rdx, rdy);

        vec sx, sy, sz;

        sx = rox - v0x;
        sy = roy - v0y;
        sz = roz - v0z;

        vec qx, qy, qz;

        vec vone  = vec(1.0);
        vec vzero = vec(0.0);

        vec a     = dot4(px, py, pz, e1x, e1y, e1z);

        // `slash dot' operator means approximate division if possible.
        vec rpa   = vone /. a;

        // q = s x e1
        qx = cross4(e1z, e1y, sy, sz);
        qy = cross4(e1x, e1z, sz, sx);
        qz = cross4(e1y, e1x, sx, sy);

        vec u, v, t;

        u = dot4(sx, sy, sz, px, py, pz) * rpa;
        v = dot4(rdx, rdy, rdz, qx, qy, qz) * rpa;
        t = dot4(e2x, e2y, e2z, qx, qy, qz) * rpa;

        vec eps;

        eps = vec(0.00001);

        vec mask0;

        mask0 = (((a * a) > eps) & ((u + v)  vzero) & (v > vzero));

        vec mask;

        //
        // final mask
        //
        mask = (mask0 & t) & (t < prevT);

        outT = sel(outT, t, mask);
        outU = sel(outU, t, mask);
        outV = sel(outV, t, mask);

        return mask;
}

MUDA compiler compiles this code into…


$ mudah isect.mu

...

__m128 cross4 (__m128 a, __m128 b, __m128 c, __m128 d)
{
    const __m128 t_vec1  =  _mm_mul_ps( a ,  b  ) ;
    const __m128 t_vec2  =  _mm_mul_ps( c ,  d  ) ;
    const __m128 t_vec3  =  _mm_sub_ps( t_vec1 ,  t_vec2  ) ;
    return t_vec3 ;
}

__m128 dot4 (__m128 ax, __m128 ay, __m128 az, __m128 bx, __m128 by, __m128 bz)
{
    const __m128 t_vec4  =  _mm_mul_ps( ax ,  bx  ) ;
    const __m128 t_vec5  =  _mm_mul_ps( ay ,  by  ) ;
    const __m128 t_vec6  =  _mm_add_ps( t_vec4 ,  t_vec5  ) ;
    const __m128 t_vec7  =  _mm_mul_ps( az ,  bz  ) ;
    const __m128 t_vec8  =  _mm_add_ps( t_vec6 ,  t_vec7  ) ;
    return t_vec8 ;
}

__m128 isect (__m128 rox, __m128 roy, __m128 roz,
__m128 rdx, __m128 rdy, __m128 rdz,
__m128 v0x, __m128 v0y, __m128 v0z,
__m128 e1x, __m128 e1y, __m128 e1z,
__m128 e2x, __m128 e2y, __m128 e2z,
__m128 * outT, __m128 * outU, __m128 * outV, __m128 prevT)
{
    __m128 px ;

    __m128 py ;

    __m128 pz ;

    const __m128 t_vec9 = cross4 (e2z, e2y, rdy, rdz) ;
    px = t_vec9 ;

    const __m128 t_vec10 = cross4 (e2x, e2z, rdz, rdx) ;
    py = t_vec10 ;

    const __m128 t_vec11 = cross4 (e2y, e2x, rdx, rdy) ;
    pz = t_vec11 ;

    __m128 sx ;

    __m128 sy ;

    __m128 sz ;

    const __m128 t_vec12  =  _mm_sub_ps( rox ,  v0x  ) ;
    sx = t_vec12 ;

    const __m128 t_vec13  =  _mm_sub_ps( roy ,  v0y  ) ;
    sy = t_vec13 ;

    const __m128 t_vec14  =  _mm_sub_ps( roz ,  v0z  ) ;
    sz = t_vec14 ;

    __m128 qx ;

    __m128 qy ;

    __m128 qz ;

    __m128 vone ;

    __m128 vzero ;

    __m128 a ;

    __m128 rpa ;

    const __m128 t_vec15 = cross4 (e1z, e1y, sy, sz) ;
    qx = t_vec15 ;

    const __m128 t_vec16 = cross4 (e1x, e1z, sz, sx) ;
    qy = t_vec16 ;

    const __m128 t_vec17 = cross4 (e1y, e1x, sx, sy) ;
    qz = t_vec17 ;

    __m128 u ;

    __m128 v ;

    __m128 t ;

    const __m128 t_vec18 = dot4 (sx, sy, sz, px, py, pz) ;
    const __m128 t_vec19  =  _mm_mul_ps( t_vec18 ,  rpa  ) ;
    u = t_vec19 ;

    const __m128 t_vec20 = dot4 (rdx, rdy, rdz, qx, qy, qz) ;
    const __m128 t_vec21  =  _mm_mul_ps( t_vec20 ,  rpa  ) ;
    v = t_vec21 ;

    const __m128 t_vec22 = dot4 (e2x, e2y, e2z, qx, qy, qz) ;
    const __m128 t_vec23  =  _mm_mul_ps( t_vec22 ,  rpa  ) ;
    t = t_vec23 ;

    __m128 eps ;

    const float t_float24 = 1.0e-5 ;
    const __m128 t_vec25 = _mm_set_ps1(  t_float24  ) ;
    const __m128 t_vec26 = t_vec25 ;
    eps = t_vec26 ;

    __m128 mask0 ;

    const __m128 t_vec27  =  _mm_mul_ps( a ,  a  ) ;
    const __m128 t_vec28  =  _mm_cmpgt_ps( t_vec27 ,  eps  ) ;
    const __m128 t_vec29  =  _mm_add_ps( u ,  v  ) ;
    const __m128 t_vec30  =  _mm_cmplt_ps( t_vec29 ,  vone  ) ;
    const __m128 t_vec31  =  _mm_and_ps( t_vec28 ,  t_vec30  ) ;
    const __m128 t_vec32  =  _mm_cmpgt_ps( u ,  vzero  ) ;
    const __m128 t_vec33  =  _mm_cmpgt_ps( v ,  vzero  ) ;
    const __m128 t_vec34  =  _mm_and_ps( t_vec32 ,  t_vec33  ) ;
    const __m128 t_vec35  =  _mm_and_ps( t_vec31 ,  t_vec34  ) ;
    mask0 = t_vec35 ;

    __m128 mask ;

    const __m128 t_vec36  =  _mm_and_ps( mask0 ,  t  ) ;
    const __m128 t_vec37  =  _mm_cmplt_ps( t ,  prevT  ) ;
    const __m128 t_vec38  =  _mm_and_ps( t_vec36 ,  t_vec37  ) ;
    mask = t_vec38 ;

    const __m128 t_vec39 = muda_sel_ps ((*outT), t, mask) ;
    (*outT) = t_vec39 ;

    const __m128 t_vec40 = muda_sel_ps ((*outU), t, mask) ;
    (*outU) = t_vec40 ;

    const __m128 t_vec41 = muda_sel_ps ((*outV), t, mask) ;
    (*outV) = t_vec41 ;

    return mask ;
}

TODOs

Here is TODOs I’m planning until the end of this year.

– Documentation!
– Support for struct
– Support for C/C++ integration(header file generation)
CodeGen for LLVM IR
– swizzling. e.g. vec dst = src.xxyy
– Support for integer and int x 4, double x 4 and double x 2 SIMD.
– Integration with AMD simulation utilities(For counting instructions)

[Ja]

とりあえず Haskell で MUDA コンパイラ書き直しました。
OCaml で初期実装したときの状態くらいにまではなっています。

ここまで時間がかかったのは、シンボルテーブルの管理でどうしても逐次処理 & 大域変数処理を伴うので、
それをどう Haskell にマッピングするかをうんうん考えていたからでした。

Hashtable, seq, unsafePerformIO とかいろいろ試したけど、結局のところ State モナドを使うことに。
Haskell 習得(主にモナド)も含めてですが、そんなんで一ヶ月半くらい費やしてしまった。

インド工科大(IIT)生なら一週間で C 言語を習得するとのことですから、
私ももっと精進しないとですね。

ついでにいうと、OCaml なら破壊的代入ですぐ解決なんですけどね。
まあ関数言語的な考え方はそのぶん習得できましたので、よしとしたいところ。

Haskell の感想

以下、 Haskell の感想

– 状態や逐次処理性を持つアルゴリズム、たとえばシンボルテーブルの管理とかがとっても書きづらい。

State モナドを使っているが、局所的に State モナドを持たせるのができなくて、
関数の上の方まで State モナドの型を”感染”させなくてはならない。
なので結局 do ばっかりのコードになる。

とはいえ、私が単に関数言語的な考え方ができていないので、
Haskell 書きづらいと感じているだけなのかもしれません。
関数型脳な人ならシンボルテーブル管理ももっとじつはうまい方法を思いつくのやも。

C などの手続き言語から入った手続き脳な私は、
どうしてもアルゴリズムを手続きで考えてしまって、
それを関数型にマッピングするという思考になってしまいます。
関数型脳で物ごとを考えるチカラがほしいです。

– コンパイルが遅い

C++ コンパイラ並みにコンパイルが遅い(GHC コンパイラの場合)…

– やっぱコンパイラ書くなら OCaml のほうが書きやすいかも。

だいたいコンパイラは逐次処理が主になるので、遅延評価がかえって邪魔になる。
しかも最近は python や ruby でも関数型的なコードが書けるので、
軽い DSL 向けコンパイラだったら python や ruby を使うのが実は一番手っ取り早いかも。

– cabal(configure みたいなの), darcs(svn みたいなの) のツールチェインで開発は魅力的

とはいえ、Haskell にはいいところもあって、OCaml よりはコードがきれいに書ける(気がする)のと、

http://haskell.org/haskellwiki/How_to_write_a_Haskell_program

にあるように、Cabal でお手軽にプロジェクトのパッケージング&ビルド環境が作れて、
darcs で個人でもお手軽にリビジョン管理できて、
QuickCheck でお手軽にランダムテストが出来て、
haddock でソースコードドキュメンテーションができる。
というよく出来た開発ツールチェインが使えるのはよい。
これだけでも魅力的で Haskell を使おうという理由になる。
# darcs はパッチが多くなると処理がもっさりするのが難点ですが…
# まあその場合は mercurial で置き換えればいいかな.

シェーダ言語も関数型で?

いままでまわりのグラフィックス野郎は C/C++ とかの手続き型言語を利用していて、
関数型言語を知っている人が皆無だったのもあり、関数型言語は私にとって新鮮でした。

何を書きたいかにもよりますが、関数型のほうがよりすっきりとコードが書ける気がします。
(基本的にアルゴリズムそのものが書けるのがよい)

今のシェーディングアルゴリズムやシェーダ言語は、
関数型にすることですっきり書けるような場面が多いのではないかという気がしています。

もちろん、関数型は手続き型にくらべて実行時のパフォーマンスに懸念がある、
という思う方が多いかと思いますが、
これは手続き型から入った人間の先入観でしかないのかもしれません。

DSL(Domain Specific Language)アプローチで適用対象に制限をもうけて、
グラフィックス向けの関数型言語なら、それ向けのコンパイラが頑張ればパフォーマンスには
問題がなくなるのではないかと思っています。

まあグラフィックス用言語に関数型がありえるとしても、
しばらくは python みたいに基本手続き型だけど、部分的に高階関数やλ式などの
関数型な言語仕様を取り込む、みたいなやり方からでしょうかね。

[1] Idea: MUDA, MUltiple Data Accelerator language for high performance computing
http://lucille.atso-net.jp/blog/?p=322

[2] Load to MUDA. SIMD code generation, domain specific language, automatic optimization, functional programming, etc.
http://lucille.atso-net.jp/blog/?p=323

[3] It’s time to move to Haskell to code MUDA?
http://lucille.atso-net.jp/blog/?p=336

MINILIGHT, a minimal global illumination renderer

MINILIGHT, a minimal global illumination renderer


http://www.hxa7241.org/minilight/minilight.html

C++ (950 lines)
OCaml (462 lines)
Ruby (499 lines)
Lua (566 lines)
Adobe Flex/ActionScript3 (642 equivalent lines (862 including extra UI))

AS でパストレかよ!というのもありますが、
もうひとつ注目したのは各言語のコード行数。

関数型言語である OCaml が一番コード行数が少ないですね。
そして実行速度は C++ > OCaml > スクリプト言語の順番とのこと。
手続き型に比べて関数型はコード量が少なくて済んで、かつ速度もそこそこ、
という利点が示されたよい例だと思います。

もっとグラフィックス野郎にも関数型言語が普及してもよいのではないかと思います。
もちろん、関数型言語でも得手不得手はありますから、
「関数型が適していそうな場所なのに、C/C++ が使われている場所があるなら、
そこは関数型に置き換えるのもいいんじゃない?」という意味です。

関数型言語の特徴は、「OCaml のすすめ」[1] に集約されています。
Ocaml, Haskell と関数型言語を一通り触ってみて、ここに書いてあることは
まったくもって当てはまるなぁと感じます。

[1] OCaml プログラミング入門
http://www.i.kyushu-u.ac.jp/~bannai/ocaml-intro/

RE: Real Time Ray-Tracing: The End of Rasterization?

http://ompf.org/forum/viewtopic.php?t=617

れしぇとふがコメントしてますね。記事にバイアスがかかっていた理由が納得。
一つの企業なのに情報の交換がなかったり(外に向けての)意見が違っているのは、大企業病らしくていいですね。
とはいえ、Intel の blog は良い意味で煽っているということでレイトレが発展することに期待。

# 本当は横目に見ているんじゃなくて、自分でも道を切り開いていくべきなのですが、
# ompf のホントにすげー奴らを見ていると、
# まだまだ自分は精進すべき段階だなぁと思います。
# だって世界中のグラフィックス研究者(主にレイトレ野郎だけど)が参加してしていて、
# フォーラムでフランクに発言しているんですよ! 世の中変わったなぁ