TextMeshPro の隠しタグ一覧

はじめに

TextMeshPro には、 ヘルプ に載っていないタグやパラメータが一部存在します。
それらを紹介していきたいと思います。

本稿の内容は TextMeshPro 3.0.6 で確認しています。

ちなみに、タグの長さの上限は 128 文字です。 <link> タグなどで長い URL を直指定した場合問題になることがあります。

一覧

<br>

改行します。 \u000a と同じです。

<nbsp>

ノーブレークスペースを入れます。 \u00a0 と同じです。

<zwsp>

ゼロ幅スペースを入れます。 \u200b と同じです。

<mark>

color および padding パラメータが使用可能です。

color は略記法 (<mark=#ff000088>) が公式ヘルプには載っていますが、明示的に指定することもできます。

padding は読んで字のごとくパディングを設定します。
<mark padding=左,右,上,下> のように指定します。負の値も指定できるようです。

<material>

マテリアルを指定します。 <material="..."> のように指定します。
実質 <font> の material パラメータでできることと同じな気がします。

</a>

パースはされますが、実装されていません。
閉じタグだけパースされます。

<gradient>

tint パラメータが使用できます。 0 か 1 を指定するようです。
<gradient tint=1> のように指定します。詳しくは試してみてください。（丸投げ）

<class>

ソースコード中に痕跡がありますが、実装されておらず、パースもされません。

<sprite>

anim パラメータが指定可能です。アニメーションします。
<sprite=1 anim=0,16,12> と指定すると、スプライト 0 番から 16 番までを 12 fps で再生します。
ランタイムでのみアニメーションするので注意してください。

<margin>

left および right パラメータが指定可能です。
それぞれ <margin-left> ・ <margin-right> と同じ働きになります。

<action>

デバッグログが表示されます。それだけです。
おそらく未実装です。表示された瞬間に何らかのコールバックを呼べるようにする予定だったのでしょうか？

<scale>

<scale=2>hoge</scale> とすると、 hoge が横に 2 倍の長さで描画されます。
なんで…？

</table> ・ <tr> ・ <th>

パースはされますが実装されていません。
</table> は閉じタグだけ認識されます。
なお、 <td> はソースコード上に痕跡がありますが実装されていません。

おわりに

<sprite anim=...> を紹介したかっただけです。
なお、ソースコードを読んだ雰囲気からすると、 TODO や懸念点が書いてあったりするため、これらが正式にサポートされるかは微妙なところです。
自己責任でお使いください。

TextMeshPro でピクセルフォントの影・縁取りをつける

はじめに

TextMeshPro は、 Unity でテキストを扱ううえでもはやなくてはならない存在です。
その表現力の高さは皆さんもご存知でしょう。

ところが、ピクセルフォントを扱うとなると数々の難がある、ということはあまり知られていないかもしれません。

影・縁取りをつけたいだけなのに

TextMeshPro の通常のフォント形式 (SDF) では、先ほど上げた画像のように影 (Underlay) や縁取り (Outline) などを自由につけることができます。

しかし、ピクセルフォントを扱う形式 (RASTER) では、インスペクタを見てもそういった項目がほとんどありません。色とテクスチャしか変更できません。
これはシェーダと描画方式が異なるためです。ピクセルフォントは TextMeshPro/Bitmap シェーダで描画されます。対して、 SDF フォントは TextMeshPro/Distance Field で描画されています。前述の影や縁取りなどは後者のシェーダでしか利用できません。

ゲームにおいて、縁取りや影は重要です。見た目のリッチさはもちろん、視認性にも寄与します。
にもかかわらず、それらの機能が使えないとなるとしんどいです。
一応、オブジェクトを複数個重ねて無理やり再現することもできなくはないですが、当然負荷も n 倍になりますし管理も面倒です。できればシェーダ側の機能でどうにかしたいところ。
また、 SDF フォントで無理やり扱う方法もなくはないですが、当然ピクセル感は失われますし、ピクセルパーフェクトな縁取りや影を付けるのは困難です。

実装要件

• ピクセルフォントを使いたい。
• それに縁取りや影を付けたい。
  • ちゃんとピクセルパーフェクトのまま縁取りや影を付けたい。
• ローカライゼーション対応の為、フォントアセットも分割しておきたい。
  • 日本語と英語への対応を仮定

なお、今回ピクセルフォントには マルモニカ を使用させていただきます。
文字サイズは x12y16 (px) です。

やったこと

フォントアセットの作成

まずは、こちらのブログを参考に、フォントアセットを作成します。

https://blog.kyubuns.dev/entry/2021/02/06/001609

フォントをインポートしたら、 Font Asset Creator を開きます。

先のブログと違う点としては、今回は滑らかな SDF ではなくピクセルの RASTER を利用するため、設定値が異なります。具体的には、

• ベースフォント
  • Sampling Point Size: Custom Size, 16
    • 必ずフォントサイズと同一にします。
  • Padding: 2
  • Packing Method: Optimum
  • Atlas Resolution: 256x256
  • Character Set: Unicode Range (Hex)
    • 詳細は後述
  • Render Mode: RASTER_HINTED
• 日本語フォント (ベースフォントと同じ項目は省略)
  • Atlas Resolution: 1024x1024
• ダイナミックフォント (〃)
  • Atlas Resolution: 1024x1024

Character Sequence には、以下を指定します。

ベースフォント:

20-7E,A0-FF,2000-200B,2010-2027,202F-2040,2044,2047-2049,20AC,2122,25A1  

日本語フォント:

A7-A8,B0-B1,B4,B6,D7,F7,2003,2010,2015-2016,2020-2021,2025-2026,2030,2032-2033,203B,2103,212B,2190-2193,21D2,21D4,2200,2202-2203,2207-2208,220B,2212,221A,221D-221E,2220,2227-222C,2234-2235,223D,2252,2260-2261,2266-2267,226A-226B,2282-2283,2286-2287,22A5,2312,25A0-25A1,25B2-25B3,25BC-25BD,25C6-25C7,25CB,25CE-25CF,25EF,2605-2606,2640,2642,266A,266D,266F,3000-3003,3005-3015,301C,3041-3093,309B-309E,30A1-30F6,30FB-30FE,4E00-4E01,4E03,4E07-4E0B,4E0D-4E0E,4E14,4E16,4E18-4E19,4E21,4E26,4E2D,4E32,4E38-4E39,4E3B-4E3C,4E45,4E4F,4E57,4E59,4E5D-4E5E,4E71,4E73,4E7E,4E80,4E86,4E88-4E89,4E8B-4E8C,4E92,4E94-4E95,4E9C,4EA1,4EA4,4EAB-4EAD,4EBA,4EC1,4ECA-4ECB,4ECF,4ED5-4ED6,4ED8-4ED9,4EDD,4EE3-4EE5,4EEE,4EF0,4EF2,4EF6,4EFB,4F01,4F0E-4F11,4F1A,4F1D,4F2F,4F34,4F38,4F3A,4F3C,4F46,4F4D-4F50,4F53,4F55,4F59,4F5C,4F73,4F75,4F7F,4F8B,4F8D,4F9B,4F9D,4FA1,4FAE-4FAF,4FB5-4FB6,4FBF,4FC2-4FC3,4FCA,4FD7,4FDD,4FE1,4FEE,4FF3,4FF5,4FF8,4FFA,5009,500B,500D,5012,5019,501F,5023-5024,502B,5039,5049,504F,505C,5065,5074-5076,507D,508D,5091,5098-5099,50AC,50B2,50B5,50B7,50BE,50C5,50CD,50CF,50D5,50DA,50E7,5100,5104,5112,511F,512A,5143-5146,5148-5149,514B,514D,5150,515A,5165,5168,516B-516D,5171,5175,5177-5178,517C,5185-5186,518A,518D,5192,5197,5199,51A0,51A5,51AC,51B6-51B7,51C4,51C6,51CD,51DD,51E1,51E6,51F6,51F8-51FA,5200,5203,5206-5208,520A,5211,5217,521D,5224-5225,5229,5230,5236-523B,5247,524A,524D,5256,525B,525D,5263-5264,526F-5270,5272,5275,5287,529B,529F-52A0,52A3,52A9-52AA,52B1,52B4,52B9,52BE,52C3,52C5,52C7,52C9,52D5,52D8-52D9,52DD,52DF,52E2,52E4,52E7,52F2,52FE,5302,5305,5316-5317,5320,5339-533B,533F,5341,5343,5347-5348,534A,5351-5354,5357-5358,535A,5360,5370-5371,5373-5375,5378,5384,5398,539A,539F,53B3,53BB,53C2,53C8,53CA-53CE,53D4,53D6-53D7,53D9,53E3-53E5,53EB-53EC,53EF-53F0,53F2-53F3,53F7-53F8,5404,5408-5409,540C-5411,541B,541F,5426,542B,5438-5439,5442,5448-544A,5468,546A,5473,547C-547D,548C,54B2,54BD,54C0-54C1,54E1,54F2,54FA,5504,5506-5507,5510,552F,5531,553E,5546,554F,5553,5584,5589,559A,559C-559D,55A9-55AB,55B6,55C5,55E3,5606,5631-5632,5668,5674,5687,56DA-56DB,56DE,56E0,56E3,56F0,56F2-56F3,56FA,56FD,570F,5712,571F,5727-5728,5730,5742,5747,574A,5751,576A,5782,578B,57A3,57CB,57CE,57DF,57F7,57F9-57FA,57FC,5800,5802,5805-5806,5815,5824,582A,5831,5834,5840-5841,584A,5851,5854,5857,585A,585E,5861,5869,587E,5883,5893,5897,589C,58A8,58B3,58BE,58C1,58C7,58CA,58CC,58EB,58EE,58F0-58F2,5909,590F,5915-5916,591A,591C,5922,5927,5929-592B,592E,5931,5947-5949,594F,5951,5954,5965,5968,596A,596E,5973-5974,597D,5982-5984,598A,5996,5999,59A5,59A8,59AC,59B9,59BB,59C9,59CB,59D3-59D4,59EB,59FB,59FF,5A01,5A18,5A20,5A2F,5A46,5A5A,5A66,5A7F,5A92,5A9B,5AC1,5AC9,5ACC,5AE1,5B22,5B50,5B54,5B57-5B58,5B5D,5B63-5B64,5B66,5B6B,5B85,5B87-5B89,5B8C,5B97-5B9D,5B9F,5BA2-5BA4,5BAE,5BB0,5BB3-5BB6,5BB9,5BBF,5BC2,5BC4,5BC6,5BCC,5BD2,5BDB,5BDD,5BDF,5BE1,5BE7,5BE9,5BEE,5BF8,5BFA,5BFE-5BFF,5C01-5C02,5C04,5C06,5C09-5C0B,5C0E-5C0F,5C11,5C1A,5C31,5C3A-5C40,5C45,5C48,5C4A-5C4B,5C55,5C5E,5C64-5C65,5C6F,5C71,5C90,5CA1,5CA9,5CAC,5CB3,5CB8,5CE0-5CE1,5CF0,5CF6,5D07,5D0E,5D16,5D29,5D50,5DDD-5DDE,5DE1,5DE3,5DE5-5DE8,5DEE,5DF1,5DFB,5DFE,5E02-5E03,5E06,5E0C,5E1D,5E25,5E2B,5E2D,5E2F-5E30,5E33,5E38,5E3D,5E45,5E55,5E63,5E72-5E74,5E78-5E79,5E7B-5E7E,5E81,5E83,5E8A,5E8F,5E95,5E97,5E9C,5EA6-5EA7,5EAB,5EAD,5EB6-5EB8,5EC3,5EC9-5ECA,5EF6-5EF7,5EFA,5F01,5F04,5F0A,5F0F-5F10,5F13-5F15,5F1F,5F25-5F27,5F31,5F35,5F37,5F3E,5F53,5F59,5F62,5F69,5F6B,5F70-5F71,5F79,5F7C,5F80-5F81,5F84-5F85,5F8B-5F8C,5F90,5F92-5F93,5F97,5FA1,5FA9-5FAA,5FAE,5FB3-5FB4,5FB9,5FC3,5FC5,5FCC-5FCD,5FD7-5FD9,5FDC,5FE0,5FEB,5FF5,6012,6016,601D,6020,6025,6027-6028,602A,604B,6050,6052,6063,6065,6068-6069,606D,606F,6075,6094,609F-60A0,60A3,60A6,60A9-60AA,60B2,60BC,60C5,60D1,60DC,60E7-60E8,60F0,60F3,6101,6109,610F,611A-611B,611F,6144,6148,614B-614C,614E,6155,6162-6163,6168,616E,6170,6176,6182,618E,61A4,61A7,61A9,61AC,61B2,61B6,61BE,61C7,61D0,61F2,61F8,6210-6212,621A,6226,622F,6234,6238,623B,623F-6240,6247,6249,624B,624D,6253,6255,6271,6276,6279,627F-6280,6284,628A,6291,6295,6297-6298,629C,629E,62AB,62B1,62B5,62B9,62BC-62BD,62C5,62C9,62CD,62D0,62D2-62D3,62D8-62D9,62DB,62DD,62E0-62E1,62EC-62ED,62F3,62F6-62F7,62FE,6301,6307,6311,6319,631F,6328,632B,632F,633F,6349,6355,6357,635C,6368,636E,637B,6383,6388,638C,6392,6398,639B,63A1-63A2,63A5,63A7-63A8,63AA,63B2,63CF-63D0,63DA-63DB,63E1,63EE,63F4,63FA,640D,642C-642D,643A,643E,6442,6458,6469,646F,6483,64A4,64AE,64B2,64C1,64CD,64E6,64EC,652F,6539,653B,653E-653F,6545,654F,6551,6557,6559,6562-6563,656C,6570,6574-6575,6577,6587,6589,658E,6591,6597,6599,659C,65A4-65A5,65AC-65AD,65B0,65B9,65BD,65C5,65CB,65CF,65D7,65E2,65E5-65E9,65EC,65FA,6606-6607,660E,6613-6614,661F-6620,6625,6627-6628,662D,662F,663C,6642,6669,666E-666F,6674,6676,6681,6687,6691,6696-6697,66A6,66AB,66AE,66B4,66C7,66D6,66DC,66F2,66F4,66F8-66F9,66FD,66FF-6700,6708-6709,670D,6715,6717,671B,671D,671F,6728,672A-672D,6731,6734,673A,673D,6749,6750-6751,675F,6761,6765,676F,6771,677E-677F,6790,6795,6797,679A,679C-679D,67A0,67A2,67AF,67B6,67C4,67D0,67D3-67D4,67F1,67F3,67F5,67FB,67FF,6803-6804,6813,6821,682A,6838-6839,683C-683D,6841,6843,6848,6851,685C,685F,6885,6897,68A8,68B0,68C4,68CB,68D2,68DA,68DF,68EE,68FA,6905,690D-690E,691C,696D,6975,6977,697C-697D,6982,69CB,69D8,69FD,6A19,6A21,6A29-6A2A,6A39,6A4B,6A5F,6B04,6B20-6B21,6B27,6B32,6B3A,6B3E,6B4C,6B53,6B62-6B63,6B66,6B69,6B6F,6B73-6B74,6B7B,6B89-6B8B,6B96,6BB4-6BB5,6BBA-6BBB,6BBF-6BC0,6BCD-6BCE,6BD2,6BD4,6BDB,6C0F,6C11,6C17,6C34,6C37-6C38,6C3E,6C41-6C42,6C4E,6C57,6C5A,6C5F-6C60,6C70,6C7A,6C7D,6C83,6C88,6C96,6C99,6CA1-6CA2,6CB3,6CB8-6CB9,6CBB-6CBC,6CBF,6CC1,6CC9-6CCA,6CCC,6CD5,6CE1-6CE3,6CE5,6CE8,6CF0,6CF3,6D0B,6D17,6D1E,6D25,6D2A,6D3B,6D3E,6D41,6D44-6D45,6D5C,6D66,6D6A,6D6E,6D74,6D77-6D78,6D88,6D99,6DAF,6DB2,6DBC,6DD1,6DE1,6DEB,6DF1,6DF7,6DFB,6E05,6E07-6E09,6E0B,6E13,6E1B,6E21,6E26,6E29,6E2C,6E2F,6E56,6E67,6E6F,6E7E-6E80,6E90,6E96,6E9D,6EB6,6EBA,6EC5,6ECB,6ED1,6EDD-6EDE,6EF4,6F01-6F02,6F06,6F0F,6F14,6F20,6F22,6F2B-6F2C,6F38,6F54,6F5C,6F5F,6F64,6F6E,6F70,6F84,6FC0-6FC1,6FC3,6FEB,6FEF,702C,706B,706F-7070,707D,7089-708A,708E,70AD,70B9-70BA,70C8,7121,7126,7136,713C,714E,7159,7167,7169,716E,718A,719F,71B1,71C3,71E5,7206,722A,7235-7236,723D,7247-7248,7259,725B,7267,7269,7272,7279,72A0,72AC,72AF,72B6,72C2,72D9,72E9,72EC-72ED,731B,731F,732B,732E,7336,733F,7344,7363,7372,7384,7387,7389,738B,73A9,73CD,73E0,73ED,73FE,7403,7406,7434,7460,7483,74A7,74B0,74BD,74E6,74F6,7518,751A,751F,7523,7528,7530-7533,7537,753A-753B,754C,754F,7551,7554,7559,755C-755D,7565,756A,7570,7573,757F,758E,7591,75AB,75B2,75BE,75C5,75C7,75D5,75D8,75DB,75E2,75E9,75F4,760D,7642,7652,7656,767A-767B,767D-767E,7684,7686-7687,76AE,76BF,76C6,76CA,76D7,76DB,76DF,76E3-76E4,76EE,76F2,76F4,76F8,76FE,7701,7709,770B-770C,771F-7720,773A,773C,7740,7761,7763,7766,77AC-77AD,77B3,77DB,77E2,77E5,77ED,77EF,77F3,7802,7814-7815,7832,7834,785D,786B-786C,7881,7891,78BA,78C1,78E8,7901,790E,793A,793C,793E,7948-7949,7956,795D-795E,7965,7968,796D,7981,7985,798D,798F,79C0-79C1,79CB,79D1-79D2,79D8,79DF,79E9,79F0,79FB,7A0B,7A0E,7A1A,7A2E,7A32,7A3C-7A3D,7A3F-7A40,7A42,7A4D,7A4F,7A6B,7A74,7A76,7A7A,7A81,7A83,7A92-7A93,7A9F,7AAE-7AAF,7ACB,7ADC,7AE0,7AE5,7AEF,7AF6,7AF9,7B11,7B1B,7B26,7B2C,7B46,7B49,7B4B,7B52,7B54,7B56,7B87,7B8B,7B97,7BA1,7BB1,7BB8,7BC0,7BC4,7BC9,7BE4,7C21,7C3F,7C4D,7C60,7C73,7C89,7C8B,7C92,7C97-7C98,7C9B,7CA7,7CBE,7CD6,7CE7,7CF8,7CFB,7CFE,7D00,7D04-7D05,7D0B,7D0D,7D14,7D19-7D1B,7D20-7D22,7D2B,7D2F-7D30,7D33,7D39-7D3A,7D42,7D44,7D4C,7D50,7D5E,7D61,7D66,7D71,7D75-7D76,7D79,7D99-7D9A,7DAD,7DB1-7DB2,7DBB,7DBF,7DCA,7DCF,7DD1-7DD2,7DDA,7DE0,7DE8-7DE9,7DEF,7DF4,7DFB,7E01,7E04,7E1B,7E26,7E2B,7E2E,7E3E,7E41,7E4A,7E54-7E55,7E6D,7E70,7F36,7F6A,7F6E,7F70,7F72,7F75,7F77,7F85,7F8A,7F8E,7F9E,7FA4,7FA8-7FA9,7FBD,7FC1,7FCC,7FD2,7FFB-7FFC,8001,8003,8005,8010,8015,8017,8033,8056,805E,8074,8077,8089,808C,8096,8098,809D,80A1-80A2,80A5,80A9-80AA,80AF,80B2,80BA,80C3,80C6,80CC,80CE,80DE,80F4,80F8,80FD,8102,8105,8107-8108,810A,811A,8131,8133,814E,8150,8155,816B,8170,8178-817A,819A,819C-819D,81A8,81B3,81C6,81D3,81E3,81E8,81EA,81ED,81F3-81F4,81FC,8208,820C,820E,8217,821E-821F,822A,822C,8236-8237,8239,8247,8266,826F,8272,8276,828B,829D,82AF,82B1,82B3,82B8,82BD,82D7,82DB,82E5-82E6,82F1,8302,830E,8328,8336,8349,8352,8358,8377,83CA,83CC,83D3,83DC,83EF,840E,843D,8449,8457,845B,846C,84B8,84C4,84CB,8511,8535,853D,8584,85A6,85AA-85AC,85CD,85E4,85E9,85FB,864E,8650,865A,865C,865E,866B,8679,868A,8695,86C7,86CD,86EE,8702,871C,878D,8840,8846,884C,8853,8857,885B,885D,8861,8863,8868,8870,8877,888B,8896,88AB,88C1-88C2,88C5,88CF,88D5,88DC,88F8,88FD-88FE,8907,8910,8912,895F,8972,897F,8981,8986-8987,898B,898F,8996,899A,89A7,89AA,89B3,89D2,89E3,89E6,8A00,8A02-8A03,8A08,8A0E,8A13,8A17-8A18,8A1F,8A2A,8A2D,8A31,8A33-8A34,8A3A,8A3C,8A50,8A54-8A55,8A5E,8A60,8A63,8A66,8A69,8A6E,8A70-8A73,8A87,8A89,8A8C-8A8D,8A93,8A95,8A98,8A9E,8AA0,8AA4,8AAC-8AAD,8AB0,8AB2,8ABF,8AC7,8ACB,8AD6,8AE6-8AE7,8AED-8AEE,8AF8,8AFE,8B00-8B01,8B04,8B0E,8B19,8B1B,8B1D,8B21,8B39,8B58,8B5C,8B66,8B70,8B72,8B77,8C37,8C46,8C4A,8C5A,8C61,8C6A,8C8C,8C9D-8C9E,8CA0-8CA2,8CA7-8CAC,8CAF,8CB4,8CB7-8CB8,8CBB-8CBC,8CBF-8CC0,8CC2-8CC4,8CC7,8CCA,8CD3,8CDB-8CDC,8CDE,8CE0,8CE2,8CE6,8CEA,8CED,8CFC,8D08,8D64,8D66,8D70,8D74,8D77,8D85,8D8A,8DA3,8DB3,8DDD,8DE1,8DEF,8DF3,8DF5,8E0A,8E0F,8E2A,8E74,8E8D,8EAB,8ECA,8ECC-8ECD,8ED2,8EDF,8EE2,8EF8,8EFD,8F03,8F09,8F1D,8F29-8F2A,8F38,8F44,8F9B,8F9E,8FA3,8FB1-8FB2,8FBA,8FBC,8FC5,8FCE,8FD1,8FD4,8FEB,8FED,8FF0,8FF7,8FFD,9000-9001,9003,9006,900F-9010,9013-9014,901A,901D,901F-9020,9023,902E,9031-9032,9038,9042,9045,9047,904A-904B,904D-904E,9053-9055,905C,9060-9061,9063,9069,906D-906E,9075,9077-9078,907A,907F,9084,90A3,90A6,90AA,90B8,90CA,90CE,90E1,90E8,90ED,90F5,90F7,90FD,914C-914E,9152,9154,9162,916A,916C,9175,9177-9178,9192,919C,91B8,91C7-91C8,91CC-91CF,91D1,91DC-91DD,91E3,920D,9234,9244,925B,9262,9271,9280,9283,9285,9298,92AD,92ED,92F3,92FC,9320,9326,932C,932E-932F,9332,934B,935B,9375,938C,9396,93AE,93E1,9418,9451,9577,9580,9589,958B,9591,9593,95A2-95A3,95A5,95B2,95C7,95D8,961C,962A,9632,963B,9644,964D,9650,965B,9662-9665,966A,9670,9673,9675-9676,9678,967A,967D,9685-9686,968A,968E-968F,9694,9699,969B-969C,96A0,96A3,96B7,96BB,96C4-96C7,96CC,96D1,96E2-96E3,96E8,96EA,96F0,96F2,96F6-96F7,96FB,9700,9707,970A,971C,9727,9732,9752,9759,975E,9762,9769,9774,97D3,97F3,97FB,97FF,9802-9803,9805-9806,9808,9810-9813,9818,982D,9830,983B-983C,984C-984E,9854-9855,9858,985E,9867,98A8,98DB,98DF,98E2,98EF,98F2,98FC-98FE,9905,990A,990C,9913,9928,9996,9999,99AC,99C4-99C6,99D0,99D2,9A0E,9A12-9A13,9A30,9A5A,9AA8,9AB8,9AC4,9AD8,9AEA,9B31,9B3C,9B42,9B45,9B54,9B5A,9BAE,9BE8,9CE5,9CF4,9D8F,9DB4,9E7F,9E93,9E97,9EA6,9EBA-9EBB,9EC4,9ED2,9ED9,9F13,9F3B,9F62,FF01,FF03-FF06,FF08-FF5E,FF61-FF9F,FFE0-FFE3,FFE5  

ダイナミックフォントは当然空にしておきます。

作り終わったら、ダイナミックフォントのアセットのインスペクタから、 Atlas Population Mode を Dynamic にし、 Multi Atlas Texture と Create Dynamic Data を有効にしておきます。

そうしたら、ベースフォントのアセットから Fallback Font Assets に日本語フォント・ダイナミックフォントのアセットを設定しておきます。

先のブログのほうではこれを動的に設定していますが、今は簡単のためいったん直接設定とします。動的設定の方法については先のブログをご参照ください。

ここまで終わったら、ためしに Text (TMP) を作ってみてテストしてみましょう。

こんな感じに表示されたら OK です。
「杖で叩く」のがミソです。これらの漢字は日本語リストに入っていないので、ダイナミックフォント経由で描画されます。

影・縁取りを描く

さて、ちゃんと描画されたら、今度は影と縁取りの描画のほうに進みましょう。

前述したように、 RASTER で作られたフォントのシェーダは TextMeshPro/Bitmap となり、これの機能が貧弱なために色を変えることぐらいしかできません。

したがって、これを改造して影を付けられるようにします。

Assets/TextMesh Pro/Shaders/TMP_Bitmap がシェーダファイルなので、これをコピーして適当な場所に配置します。

そうしたら、以下のように編集します。

Shader "TextMeshPro/Bitmap_Shadow" {        // この名前を変える(重複しないように)  
  
Properties {  
    _MainTex        ("Font Atlas", 2D) = "white" {}  
    _FaceTex        ("Font Texture", 2D) = "white" {}  
    [HDR]_FaceColor ("Text Color", Color) = (1,1,1,1)  
    [HDR]_ShadowColor ("Shadow Color", Color) = (0.5,0.5,0.5,1)     // 追加  
    _ShadowFlag     ("Shadow Flag", int) = 255      // 追加  
  
    _VertexOffsetX  ("Vertex OffsetX", float) = 0  
    _VertexOffsetY  ("Vertex OffsetY", float) = 0  
    _MaskSoftnessX  ("Mask SoftnessX", float) = 0  
    _MaskSoftnessY  ("Mask SoftnessY", float) = 0  
  
    _ClipRect("Clip Rect", vector) = (-32767, -32767, 32767, 32767)  
  
    _StencilComp("Stencil Comparison", Float) = 8  
    _Stencil("Stencil ID", Float) = 0  
    _StencilOp("Stencil Operation", Float) = 0  
    _StencilWriteMask("Stencil Write Mask", Float) = 255  
    _StencilReadMask("Stencil Read Mask", Float) = 255  
  
    _CullMode("Cull Mode", Float) = 0  
    _ColorMask("Color Mask", Float) = 15  
}  
  
SubShader{  
  
    Tags { "Queue" = "Transparent" "IgnoreProjector" = "True" "RenderType" = "Transparent" }  
  
    Stencil  
    {  
        Ref[_Stencil]  
        Comp[_StencilComp]  
        Pass[_StencilOp]  
        ReadMask[_StencilReadMask]  
        WriteMask[_StencilWriteMask]  
    }  
  
  
    Lighting Off  
    Cull [_CullMode]  
    ZTest [unity_GUIZTestMode]  
    ZWrite Off  
    Fog { Mode Off }  
    Blend SrcAlpha OneMinusSrcAlpha  
    ColorMask[_ColorMask]  
  
    Pass {  
        CGPROGRAM  
        #pragma vertex vert  
        #pragma fragment frag  
  
        #pragma multi_compile __ UNITY_UI_CLIP_RECT  
        #pragma multi_compile __ UNITY_UI_ALPHACLIP  
  
  
        #include "UnityCG.cginc"  
  
        struct appdata_t {  
            float4 vertex       : POSITION;  
            fixed4 color        : COLOR;  
            float2 texcoord0    : TEXCOORD0;  
            float2 texcoord1    : TEXCOORD1;  
        };  
  
        struct v2f {  
            float4  vertex      : SV_POSITION;  
            fixed4  color       : COLOR;  
            fixed4  color2      : COLOR1;       // 追加  
            float2  texcoord0   : TEXCOORD0;  
            float2  texcoord1   : TEXCOORD1;  
            float4  mask        : TEXCOORD2;  
            float2  texcoord3   : TEXCOORD3;    // 追加  
        };  
  
        uniform sampler2D   _MainTex;  
        uniform float4      _MainTex_TexelSize;     // 追加  
        uniform sampler2D   _FaceTex;  
        uniform float4      _FaceTex_ST;  
        uniform fixed4      _FaceColor;  
        uniform fixed4      _ShadowColor;       // 追加  
        uniform int         _ShadowFlag;        // 追加  
  
        uniform float       _VertexOffsetX;  
        uniform float       _VertexOffsetY;  
        uniform float4      _ClipRect;  
        uniform float       _MaskSoftnessX;  
        uniform float       _MaskSoftnessY;  
  
        float2 UnpackUV(float uv)  
        {  
            float2 output;  
            output.x = floor(uv / 4096);  
            output.y = uv - 4096 * output.x;  
  
            return output * 0.001953125;  
        }  
  
        v2f vert (appdata_t v)  
        {  
            float4 vert = v.vertex;  
            vert.x += _VertexOffsetX;  
            vert.y += _VertexOffsetY;  
  
            vert.xy += (vert.w * 0.5) / _ScreenParams.xy;  
  
            float4 vPosition = UnityPixelSnap(UnityObjectToClipPos(vert));  
  
            fixed4 faceColor = v.color;  
            faceColor *= _FaceColor;  
            fixed4 shadowColor = _ShadowColor;          // 追加  
  
            v2f OUT;  
            OUT.vertex = vPosition;  
            OUT.color = faceColor;  
            OUT.texcoord0 = v.texcoord0;  
            OUT.texcoord1 = TRANSFORM_TEX(UnpackUV(v.texcoord1), _FaceTex);  
            float2 pixelSize = vPosition.w;  
            pixelSize /= abs(float2(_ScreenParams.x * UNITY_MATRIX_P[0][0], _ScreenParams.y * UNITY_MATRIX_P[1][1]));  
  
            // Clamp _ClipRect to 16bit.  
            float4 clampedRect = clamp(_ClipRect, -2e10, 2e10);  
            OUT.mask = float4(vert.xy * 2 - clampedRect.xy - clampedRect.zw, 0.25 / (0.25 * half2(_MaskSoftnessX, _MaskSoftnessY) + pixelSize.xy));  
  
            OUT.texcoord3 = _MainTex_TexelSize.xy;      // 追加  
            OUT.color2 = shadowColor;                   // 追加  
  
            return OUT;  
        }  
  
        fixed4 frag (v2f IN) : SV_Target  
        {  
            fixed4 color = tex2D(_MainTex, IN.texcoord0);  
            color = fixed4 (tex2D(_FaceTex, IN.texcoord1).rgb * IN.color.rgb, IN.color.a * color.a);  
  
  
            // ここから追加  
            fixed4 color2 = tex2D(_MainTex, IN.texcoord0 + float2(+IN.texcoord3.x, 0));  
            color2 = fixed4 (tex2D(_FaceTex, IN.texcoord1 + float2(+IN.texcoord3.x, 0)).rgb * IN.color2.rgb, IN.color2.a * color2.a * ((_ShadowFlag >> 0) & 1));  
  
            fixed4 color3 = tex2D(_MainTex, IN.texcoord0 + float2(-IN.texcoord3.x, 0));  
            color3 = fixed4 (tex2D(_FaceTex, IN.texcoord1 + float2(-IN.texcoord3.x, 0)).rgb * IN.color2.rgb, IN.color2.a * color3.a * ((_ShadowFlag >> 1) & 1));  
  
            fixed4 color4 = tex2D(_MainTex, IN.texcoord0 + float2(0, +IN.texcoord3.y));  
            color4 = fixed4 (tex2D(_FaceTex, IN.texcoord1 + float2(0, +IN.texcoord3.y)).rgb * IN.color2.rgb, IN.color2.a * color4.a * ((_ShadowFlag >> 2) & 1));  
  
            fixed4 color5 = tex2D(_MainTex, IN.texcoord0 + float2(0, -IN.texcoord3.y));  
            color5 = fixed4 (tex2D(_FaceTex, IN.texcoord1 + float2(0, -IN.texcoord3.y)).rgb * IN.color2.rgb, IN.color2.a * color5.a * ((_ShadowFlag >> 3) & 1));  
  
  
            fixed4 color6 = tex2D(_MainTex, IN.texcoord0 + float2(+IN.texcoord3.x, +IN.texcoord3.y));  
            color6 = fixed4 (tex2D(_FaceTex, IN.texcoord1 + float2(+IN.texcoord3.x, +IN.texcoord3.y)).rgb * IN.color2.rgb, IN.color2.a * color6.a * ((_ShadowFlag >> 4) & 1));  
  
            fixed4 color7 = tex2D(_MainTex, IN.texcoord0 + float2(-IN.texcoord3.x, -IN.texcoord3.y));  
            color7 = fixed4 (tex2D(_FaceTex, IN.texcoord1 + float2(-IN.texcoord3.x, -IN.texcoord3.y)).rgb * IN.color2.rgb, IN.color2.a * color7.a * ((_ShadowFlag >> 5) & 1));  
  
            fixed4 color8 = tex2D(_MainTex, IN.texcoord0 + float2(-IN.texcoord3.x, +IN.texcoord3.y));  
            color8 = fixed4 (tex2D(_FaceTex, IN.texcoord1 + float2(-IN.texcoord3.x, +IN.texcoord3.y)).rgb * IN.color2.rgb, IN.color2.a * color8.a * ((_ShadowFlag >> 6) & 1));  
  
            fixed4 color9 = tex2D(_MainTex, IN.texcoord0 + float2(+IN.texcoord3.x, -IN.texcoord3.y));  
            color9 = fixed4 (tex2D(_FaceTex, IN.texcoord1 + float2(+IN.texcoord3.x, -IN.texcoord3.y)).rgb * IN.color2.rgb, IN.color2.a * color9.a * ((_ShadowFlag >> 7) & 1));  
  
  
            color2 = max(max(max(max(max(max(max(color2, color3), color4), color5), color6), color7), color8), color9);  
            color = fixed4(color.rgb * color.a + color2.rgb * color2.a * (1 - color.a), color.a + color2.a * (1 - color.a));  
  
            // ここまで  
  
            // Alternative implementation to UnityGet2DClipping with support for softness.  
            #if UNITY_UI_CLIP_RECT  
                half2 m = saturate((_ClipRect.zw - _ClipRect.xy - abs(IN.mask.xy)) * IN.mask.zw);  
                color *= m.x * m.y;  
            #endif  
  
            #if UNITY_UI_ALPHACLIP  
                clip(color.a - 0.001);  
            #endif  
  
            return color;  
        }  
        ENDCG  
    }  
}  
  
    // 下の行を消す (消さないと標準 GUI が出て ShadowColor が設定できない)  
    // CustomEditor "TMPro.EditorUtilities.TMP_BitmapShaderGUI"  
}  
  

編集が終わったら、フォントアセット内のマテリアルを開き、ためしに Shader を TextMeshPro/Bitmap_Shadow に変更してみましょう。

うまくいけば、画面上のアルファベットがグレーの枠で囲われるはずです。
日本語がそのままなのはいったん目をつぶりましょう。あとで直します。

このままだとすべての文字が囲われてしまうので、いったんマテリアルを作ります。

フォントアセット内のマテリアルを開き、 Inspector の上のほう (Shader を変更できるあたり) で右クリック→ Create Material Preset を選択すると、マテリアルが生えます。
生えたマテリアルには x12y16pxMaruMonica_Outline などと命名しましょう。この時、必ずフォントアセット名から始まるように (ここでは x12y16pxMaruMonica_ から始まる名前に) してください。
同じ手順でもう一個マテリアルを生やし、今度は x12y16pxMaruMonica_Shadow と命名します。
今作ったマテリアルを Inspector で開いて、 Shadow Flag を 70 に設定します。どうして 70 なのかというと、 70=64+4+2 で、影の位置のフラグになっています。詳しくはシェーダを参照してください。
そして、今度はフォントアセット内のマテリアルの Shadow Flag を 0 にします。こうしないと全部の文字が囲われてしまいますからね。
また、それぞれのフォントアセット内のマテリアルを開いて、 Shader を TextMeshPro/Bitmap_Shadow に変更しておきます。

以上の手順がうまくいっていれば、 TextMeshPro - Text (UI) の Material Preset から、今作ったマテリアルが選択できるようになっているはずです。
ここを切り替えれば通常/囲み/影を切り替えることができます。

お疲れさまでした！
…とはいきません。

英語フォントにしか反映されない問題

上の画像を見ればわかる通り、英語の部分にしか書式が反映されておらず、フォールバック先の日本語・ダイナミックフォント (杖) はそのままになっています。

どうしてかというと、通常 SDF フォントであれば、フォールバックしてもマテリアルのパラメータが引き継がれるのですが、 RASTER フォントではなぜかその機能がありません。くそ！

例えば、 Noto Sans JP で同様の手順で SDF フォントを三つ (ベース/日本語/ダイナミックで) 作ってベースフォントの Outline を変更すると、ちゃんと全部にアウトラインが設定されます。

本来こうあってほしいのですが、なーぜーかー RASTER ではこの機能がありません。

これをどう解決したかというと、まずソースコードを読みます。基本ですね。
TMP_MaterialManager.cs を読むと、以下の分岐が目につきます。

// Create new material from the source material and copy properties if using distance field shaders.  
Material fallbackMaterial;  
if (sourceMaterial.HasProperty(ShaderUtilities.ID_GradientScale) && targetMaterial.HasProperty(ShaderUtilities.ID_GradientScale))  

if using distance field shaders. 小癪な……
ということで、 SDF フォントかどうかを判定し、そうであれば正しくマテリアルパラメータが引き継がれるコードになっています。
これを逆手にとって、この分岐を通してさえしまえば…… シェーダが GradientScale プロパティを持ってさえいれば、正しく引継ぎが行われるはずです。

というわけで、先ほど作ったシェーダに以下のプロパティを生やします。

// 前略  
[HDR]_ShadowColor ("Shadow Color", Color) = (0.5,0.5,0.5,1)  
_ShadowFlag     ("Shadow Flag", int) = 255  
  
// ここから追加: TMP を騙せ  
_GradientScale  ("Gradient Scale", float) = 5.0       
_TextureWidth       ("Texture Width", float) = 512  
_TextureHeight      ("Texture Height", float) = 512  
_WeightNormal       ("Weight Normal", float) = 0  
_WeightBold         ("Weight Bold", float) = 0.5  
// ここまで追加  
  
_VertexOffsetX  ("Vertex OffsetX", float) = 0  
_VertexOffsetY  ("Vertex OffsetY", float) = 0  
// 後略  

_GradientScale 以外にも増えているのはほかの場所でも辻褄を合わせるためです。消すとエラーを吐くのでやってみてください。

直したところで、試してみましょう：

やりましたね！

おわりに

調査に 2 日かかりました。同じ轍を踏んで頭を抱えている人の助けになれば幸いです。

なお、シェーダ初心者のため、だいぶガバいシェーダーコードを書いているのは承知しております。マルチパスするよりはましな…はず…
こうしたほうがいいよ、などあればぜひご連絡ください。

最後におまけで、編集が完了した状態のファイルを gist に置いておきます。独自性もないので、私が変更した個所は CC0 ライセンスです。

ピクセルフォント囲み/影シェーダ

まもけん2 攻略シート

下に同じく、スプレッドシートのままだと検索に引っ掛からない可能性に思い至ったので、リンクを貼っておきます。

魔物娘と不思議な冒険2 (steam)

シートへのリンク

Cuisineer 攻略シート

スプレッドシートのままだと Google 検索に引っ掛からない可能性に思い至ったので、ここにリンクを貼っておきます。

RTA のチャートを考えたいなどあればご利用ください。

Cuisineer (steam)

シートへのリンク

Lerp 手法 5 選

はじめに

Lerp というのは 線形補間 を行う関数で、直近の 2 点間を直線で結んだとき、位置 (の比率) t にある点の値を返します。

C# では System.Numerics.Vector3.Lerp や .NET 8 から実装される Double.Lerp があります。
これらよりは Unity の Mathf.Lerp のほうがなじみ深いかもしれません。

これを行う手法は複数あり、それぞれ特徴があります。それらの紹介と比較をしていきたいと思います。

5 つの Lerp 手法

Lerp のシグネチャは以下であるものとします。

double Lerp(double x, double y, double t) => ...  

1. x + (y - x) * t

多くの方がこれを一番最初に思いつくかと思います。
min + (max - min) * t の形式で見かけた・書いた方は結構いらっしゃるのではないでしょうか。

2. (1 - t) * x + t * y

この式には「 t == 1.0 のときに必ず y と等しくなる」という特徴があります。

これは当たり前のように思えるかもしれませんが、実は 1 番目の式 x + (y - x) * t の場合、浮動小数点数の計算誤差によって t == 1.0 でも結果が y にならない場合があります (具体例は後述します) 。
そのような挙動を避けたい場合に採用されるのがこれです。

3. fma(t, y - x, x)

急に関数が出てきました。
fma() は Fused Multiply Add といって、 fma(a, b, c) ~= a * b + c となるような関数です。
~= になっているのがポイントで、 fma は計算中に丸めを 1 回だけ行います。 a * b + c は * と + で 2 回丸めが発生しますので、それよりも精度の高い計算を行うことができます。
C# では Math.FusedMultiplyAdd() から利用することができます。

ちなみに、浮動小数点数演算を最適化してくれるサイト Herbie に 1 つめの式を入れるとこれが返ってきます。

4. fma(t, y, (1 - t) * x)

これは 2 つ目の式に対して fma を適用した式です。

5. fma(t, y, fma(-t, x, x))

最後に、 4 つ目の式を変形してさらに fma を適用した式です。
((1 - t) * x → x - t * x → -t * x + x より)
ついに fma だけになりました。

手法の比較

極端な値を与えた場合の結果 (t == 0.5)

x と y にそれぞれ極端な値 ( ±∞ とか NaN とか) を指定した場合で t == 0.5 のとき、直感的には下表のようになってほしい気持ちがあります。
なお、 max は表現可能な最大値の double.MaxValue ~= 1.8e+308 、 eps は絶対値が表現可能な最小値の double.Epsilon ~= 4.9e-324 を指します。

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN	NaN
-max	-∞	-max	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	0	+∞	NaN
-1	-∞	-9e+307	-1	-0.5	-0.5	-0.5	-0.5	0	+9e+307	+∞	NaN
-eps	-∞	-9e+307	-0.5	-eps	-0	-0	0	+0.5	+9e+307	+∞	NaN
-0	-∞	-9e+307	-0.5	-0	-0	-0	0	+0.5	+9e+307	+∞	NaN
0	-∞	-9e+307	-0.5	-0	0	0	0	+0.5	+9e+307	+∞	NaN
+eps	-∞	-9e+307	-0.5	0	0	0	+eps	+0.5	+9e+307	+∞	NaN
+1	-∞	-9e+307	0	+0.5	+0.5	+0.5	+0.5	+1	+9e+307	+∞	NaN
+max	-∞	0	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+max	+∞	NaN
+∞	NaN	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

出力結果を見比べてみましょう。
上の理想の表と違う部分は太字にしてあります。

1. x + (y - x) * t

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN	NaN
-max	NaN	-max	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	-∞	NaN	NaN
-1	NaN	-9e+307	-1	-0.5	-0.5	-0.5	-0.5	0	+9e+307	NaN	NaN
-eps	NaN	-9e+307	-0.5	-eps	-0	0	0	+0.5	+9e+307	NaN	NaN
-0	NaN	-9e+307	-0.5	-eps	0	0	+eps	+0.5	+9e+307	NaN	NaN
0	NaN	-9e+307	-0.5	-eps	0	0	+eps	+0.5	+9e+307	NaN	NaN
+eps	NaN	-9e+307	-0.5	0	0	0	+eps	+0.5	+9e+307	NaN	NaN
+1	NaN	-9e+307	0	+0.5	+0.5	+0.5	+0.5	+1	+9e+307	NaN	NaN
+max	NaN	+∞	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+max	NaN	NaN
+∞	NaN	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

これを見ると、 x = ±∞ の場合に NaN になってしまっていることがわかります。
なぜかというと、例えば x = ∞ の場合は、

∞ + (y - ∞) * t  
∞ + (-∞) * t  
∞ + (-∞)  
NaN  

となってしまうためです。

加えて、 (-max, max) の場合にオーバーフローして ±∞ になってしまっています。

また、微妙な問題ですが、 (x, y) == (-0, -0) のときに +0 を返しています。

2. (1 - t) * x + t * y

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN	NaN
-max	-∞	-max	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	0	+∞	NaN
-1	-∞	-9e+307	-1	-0.5	-0.5	-0.5	-0.5	0	+9e+307	+∞	NaN
-eps	-∞	-9e+307	-0.5	-0	-0	0	0	+0.5	+9e+307	+∞	NaN
-0	-∞	-9e+307	-0.5	-0	-0	0	0	+0.5	+9e+307	+∞	NaN
0	-∞	-9e+307	-0.5	0	0	0	0	+0.5	+9e+307	+∞	NaN
+eps	-∞	-9e+307	-0.5	0	0	0	0	+0.5	+9e+307	+∞	NaN
+1	-∞	-9e+307	0	+0.5	+0.5	+0.5	+0.5	+1	+9e+307	+∞	NaN
+max	-∞	0	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+max	+∞	NaN
+∞	NaN	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

これは 1 つめの式とは異なり、概ね直感的な結果を返しています。

これも微妙な問題ですが、 Epsilon 同士のときに 0 を返しています。 x, y の項それぞれで 0.5 倍するため、アンダーフローして 0 になってしまったためです。

3. fma(t, y - x, x)

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN	NaN
-max	NaN	-max	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	-∞	NaN	NaN
-1	NaN	-9e+307	-1	-0.5	-0.5	-0.5	-0.5	0	+9e+307	NaN	NaN
-eps	NaN	-9e+307	-0.5	-eps	-0	-0	0	+0.5	+9e+307	NaN	NaN
-0	NaN	-9e+307	-0.5	-0	0	0	0	+0.5	+9e+307	NaN	NaN
0	NaN	-9e+307	-0.5	-0	0	0	0	+0.5	+9e+307	NaN	NaN
+eps	NaN	-9e+307	-0.5	0	0	0	+eps	+0.5	+9e+307	NaN	NaN
+1	NaN	-9e+307	0	+0.5	+0.5	+0.5	+0.5	+1	+9e+307	NaN	NaN
+max	NaN	+∞	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+max	NaN	NaN
+∞	NaN	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

1 つめの式と同様、 NaN と ±∞ 、 -0 の問題があります。

4. fma(t, y, (1 - t) * x)

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN	NaN
-max	-∞	-max	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	0	+∞	NaN
-1	-∞	-9e+307	-1	-0.5	-0.5	-0.5	-0.5	0	+9e+307	+∞	NaN
-eps	-∞	-9e+307	-0.5	-0	-0	-0	-0	+0.5	+9e+307	+∞	NaN
-0	-∞	-9e+307	-0.5	-0	-0	0	0	+0.5	+9e+307	+∞	NaN
0	-∞	-9e+307	-0.5	0	0	0	0	+0.5	+9e+307	+∞	NaN
+eps	-∞	-9e+307	-0.5	0	0	0	0	+0.5	+9e+307	+∞	NaN
+1	-∞	-9e+307	0	+0.5	+0.5	+0.5	+0.5	+1	+9e+307	+∞	NaN
+max	-∞	0	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+max	+∞	NaN
+∞	NaN	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

こちらは概ね 2 つめの式と同じ雰囲気で、概ね直感的です。

5. fma(t, y, fma(-t, x, x))

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN	NaN
-max	NaN	-max	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	-9e+307	0	NaN	NaN
-1	NaN	-9e+307	-1	-0.5	-0.5	-0.5	-0.5	0	+9e+307	NaN	NaN
-eps	NaN	-9e+307	-0.5	-0	-0	-0	-0	+0.5	+9e+307	NaN	NaN
-0	NaN	-9e+307	-0.5	-0	0	0	0	+0.5	+9e+307	NaN	NaN
0	NaN	-9e+307	-0.5	0	0	0	0	+0.5	+9e+307	NaN	NaN
+eps	NaN	-9e+307	-0.5	0	0	0	0	+0.5	+9e+307	NaN	NaN
+1	NaN	-9e+307	0	+0.5	+0.5	+0.5	+0.5	+1	+9e+307	NaN	NaN
+max	NaN	0	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+9e+307	+max	NaN	NaN
+∞	NaN	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

こちらも ±∞ まわりで NaN 問題が見られます。

かわりに (-max, max) での ±∞ になっていた問題は解消しています。

極端な値を与えた場合の結果 (t == 0.0 および t == 1.0)

さて、 t == 0.0 および t == 1.0 の場合、 NaN が絡まない限り x または y がそのまま返ることが期待されます。
この場合はどうでしょうか？

t == 0.0 の場合 (上表) と t == 1.0 の場合 (下表) の理想値を示します。

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-max	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-1	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-eps	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-0	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
0	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
+eps	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
+1	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
+max	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
+∞	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN
-max	-max	-max	-max	-max	-max	-max	-max	-max	-max	-max	NaN
-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	NaN
-eps	-eps	-eps	-eps	-eps	-eps	-eps	-eps	-eps	-eps	-eps	NaN
-0	-0	-0	-0	-0	-0	-0	-0	-0	-0	-0	NaN
0	0	0	0	0	0	0	0	0	0	0	NaN
+eps	+eps	+eps	+eps	+eps	+eps	+eps	+eps	+eps	+eps	+eps	NaN
+1	+1	+1	+1	+1	+1	+1	+1	+1	+1	+1	NaN
+max	+max	+max	+max	+max	+max	+max	+max	+max	+max	+max	NaN
+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

1. x + (y - x) * t

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
-max	NaN	-max	-1	-eps	-0	0	+eps	+1	NaN	NaN	NaN
-1	NaN	-max	-1	-eps	-0	0	+eps	+1	+max	NaN	NaN
-eps	NaN	-max	-1	-eps	-0	0	+eps	+1	+max	NaN	NaN
-0	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
0	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+eps	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+1	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+max	NaN	NaN	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+∞	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN	NaN
-max	NaN	-max	-max	-max	-max	-max	-max	-max	-∞	NaN	NaN
-1	NaN	0	-1	-1	-1	-1	-1	-1	0	NaN	NaN
-eps	NaN	0	0	-eps	-eps	-eps	-eps	0	0	NaN	NaN
-0	NaN	0	0	0	0	0	0	0	0	NaN	NaN
0	NaN	0	0	0	0	0	0	0	0	NaN	NaN
+eps	NaN	0	0	+eps	+eps	+eps	+eps	0	0	NaN	NaN
+1	NaN	0	+1	+1	+1	+1	+1	+1	0	NaN	NaN
+max	NaN	+∞	+max	+max	+max	+max	+max	+max	+max	NaN	NaN
+∞	NaN	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

t == 0.0 の場合、 ±∞ が絡むときと (-max, +max) のときに NaN になってしまっています。

また t == 1.0 の場合でも、 x == ±∞ のときに NaN になってしまっています。
加えて、 max が絡む場合に、上下限に含まれていない 0 を返している箇所があります。おそらく桁落ちで情報量を失ってしまったものと思います。

2. (1 - t) * x + t * y

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
-max	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-1	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-eps	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-0	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
0	-∞	-max	-1	-eps	0	0	+eps	+1	+max	+∞	NaN
+eps	-∞	-max	-1	-eps	0	0	+eps	+1	+max	+∞	NaN
+1	-∞	-max	-1	-eps	0	0	+eps	+1	+max	+∞	NaN
+max	-∞	-max	-1	-eps	0	0	+eps	+1	+max	+∞	NaN
+∞	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN	NaN
-max	NaN	-max	-max	-max	-max	-max	-max	-max	-max	NaN	NaN
-1	NaN	-1	-1	-1	-1	-1	-1	-1	-1	NaN	NaN
-eps	NaN	-eps	-eps	-eps	-eps	-eps	-eps	-eps	-eps	NaN	NaN
-0	NaN	-0	-0	-0	-0	0	0	0	0	NaN	NaN
0	NaN	0	0	0	0	0	0	0	0	NaN	NaN
+eps	NaN	+eps	+eps	+eps	+eps	+eps	+eps	+eps	+eps	NaN	NaN
+1	NaN	+1	+1	+1	+1	+1	+1	+1	+1	NaN	NaN
+max	NaN	+max	+max	+max	+max	+max	+max	+max	+max	NaN	NaN
+∞	NaN	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

1 つめの式より NaN の量が減っており、概ね x, y それぞれの値を返せています。
それでも NaN 以外から NaN になる行・列があります。
0 * ±∞ が NaN になってしまうためです。

また、 (x, y) == (-0, -0) の場合にきちんと (?) -0 を返しています。

3. fma(t, y - x, x)

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
-max	NaN	-max	-1	-eps	-0	0	+eps	+1	NaN	NaN	NaN
-1	NaN	-max	-1	-eps	-0	0	+eps	+1	+max	NaN	NaN
-eps	NaN	-max	-1	-eps	-0	0	+eps	+1	+max	NaN	NaN
-0	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
0	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+eps	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+1	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+max	NaN	NaN	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+∞	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN	NaN
-max	NaN	-max	-max	-max	-max	-max	-max	-max	-∞	NaN	NaN
-1	NaN	0	-1	-1	-1	-1	-1	-1	0	NaN	NaN
-eps	NaN	0	0	-eps	-eps	-eps	-eps	0	0	NaN	NaN
-0	NaN	0	0	0	0	0	0	0	0	NaN	NaN
0	NaN	0	0	0	0	0	0	0	0	NaN	NaN
+eps	NaN	0	0	+eps	+eps	+eps	+eps	0	0	NaN	NaN
+1	NaN	0	+1	+1	+1	+1	+1	+1	0	NaN	NaN
+max	NaN	+∞	+max	+max	+max	+max	+max	+max	+max	NaN	NaN
+∞	NaN	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

1 つめの式 x + (y - x) * t と同じ感じで、大きな値の部分で NaN が返りがちです。

4. fma(t, y, (1 - t) * x)

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
-max	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-1	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-eps	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-0	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
0	-∞	-max	-1	-eps	0	0	+eps	+1	+max	+∞	NaN
+eps	-∞	-max	-1	-eps	0	0	+eps	+1	+max	+∞	NaN
+1	-∞	-max	-1	-eps	0	0	+eps	+1	+max	+∞	NaN
+max	-∞	-max	-1	-eps	0	0	+eps	+1	+max	+∞	NaN
+∞	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN	NaN
-max	NaN	-max	-max	-max	-max	-max	-max	-max	-max	NaN	NaN
-1	NaN	-1	-1	-1	-1	-1	-1	-1	-1	NaN	NaN
-eps	NaN	-eps	-eps	-eps	-eps	-eps	-eps	-eps	-eps	NaN	NaN
-0	NaN	-0	-0	-0	-0	0	0	0	0	NaN	NaN
0	NaN	0	0	0	0	0	0	0	0	NaN	NaN
+eps	NaN	+eps	+eps	+eps	+eps	+eps	+eps	+eps	+eps	NaN	NaN
+1	NaN	+1	+1	+1	+1	+1	+1	+1	+1	NaN	NaN
+max	NaN	+max	+max	+max	+max	+max	+max	+max	+max	NaN	NaN
+∞	NaN	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

2 つめの式 (1 - t) * x + t * y と同じ感じで、 ±∞ 以外は直感的な値を返せています。

5. fma(t, y, fma(-t, x, x))

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
-max	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
-1	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
-eps	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
-0	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
0	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+eps	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+1	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+max	NaN	-max	-1	-eps	0	0	+eps	+1	+max	NaN	NaN
+∞	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

y＼x	-∞	-max	-1	-eps	-0	0	+eps	+1	+max	+∞	NaN
-∞	NaN	-∞	-∞	-∞	-∞	-∞	-∞	-∞	-∞	NaN	NaN
-max	NaN	-max	-max	-max	-max	-max	-max	-max	-max	NaN	NaN
-1	NaN	-1	-1	-1	-1	-1	-1	-1	-1	NaN	NaN
-eps	NaN	-eps	-eps	-eps	-eps	-eps	-eps	-eps	-eps	NaN	NaN
-0	NaN	0	0	0	0	0	0	0	0	NaN	NaN
0	NaN	0	0	0	0	0	0	0	0	NaN	NaN
+eps	NaN	+eps	+eps	+eps	+eps	+eps	+eps	+eps	+eps	NaN	NaN
+1	NaN	+1	+1	+1	+1	+1	+1	+1	+1	NaN	NaN
+max	NaN	+max	+max	+max	+max	+max	+max	+max	+max	NaN	NaN
+∞	NaN	+∞	+∞	+∞	+∞	+∞	+∞	+∞	+∞	NaN	NaN
NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

やはり 1. や 3. の式と似た挙動ですが、 ±∞ が絡まないあたりについては正しい値を返せています。

まとめると、 t == 0.0 および t == 1.0 の場合には、 2. と 4. の式が直感に近い挙動をすることがわかりました。
ただ、それでも完璧に x または y の値を返すわけではない (NaN を含まない状態から NaN になりうる) 点に注意が必要です。

そんな値を与えた状態で Lerp しようとするな、はおっしゃる通りですが……

精度

前述したように、 1 つめの式 x + (y - x) * t では、 t == 1.0 でも結果が y に等しくならないことがあります。具体的に見てみましょう。

x = 456789, y = 0.1, t = 1.0 のとき、それぞれ以下の値が返ります。

1. 0.09999999997671694
2. 0.1
3. 0.09999999997671694
4. 0.1
5. 0.1

1 つめの式 x + (y - x) * t と 3 つめの式 fma(t, y - x, x) では問題が出てしまっています。
x, y を乱数にして試したところ、だいたい 16 % 程度にこの現象がみられるようです。

ほかにも結果が変わる例があったので紹介します。
最初の行が (x, y, t) 、続く行が n 番目の式、最後の行 (0. がない行) は整数化して BigInteger で計算した参考値です。

0.43840423350536384, 0.35411112619051244, 0.21474226926519646  
        0.42030294035715793  
        0.4203029403571579  
        0.42030294035715793  
        0.4203029403571579  
        0.4203029403571579  
          4203029403571579193038747032939560  
  
0.2530206785167861, 0.7429566475206353, 0.42010793232833676  
        0.45884666542827324  
        0.4588466654282733  
        0.4588466654282733  
        0.4588466654282733  
        0.4588466654282733  
          458846665428273276122038673856592  
  
0.42627189060522985, 0.28161146670562975, 0.5186179477132566  
        0.3512483984470895  
        0.35124839844708955  
        0.3512483984470895  
        0.3512483984470895  
        0.3512483984470895  
          351248398447089509911817791314340  
  
0.48581429649946484, 0.05085961859894306, 0.871859806353322  
        0.10659479525264437  
        0.10659479525264434  
        0.10659479525264436  
        0.10659479525264436  
        0.10659479525264436  
          10659479525264437688525716364684  

乱数で調べた限りでは、 4 番目と 5 番目の式は常に一致するようです。

BigInteger での計算値が正しいと仮定し、それとの差を比較してみたところ、誤差が小さい順に Lerp3 < Lerp1 < Lerp4 == Lerp5 < Lerp2 という結果が得られました。

2 番目の式 (1 - t) * x + t * y の精度が良いと聞いていたので、 3 番目の式 fma(t, y - x, x) が一番良いのが意外でした。
そういえば Herbie は 3 つめの式を返していましたね…　実際に精度が良いことが確かめられた形になります。

速度

そのまま測ったところどれも測定不能レベルで速かったので、 Random.Shared.NextDouble() を各引数に渡して若干遅らせています。

BenchmarkDotNet=v0.13.5, OS=Windows 11 (10.0.22621.1555/22H2/2022Update/SunValley2)  
12th Gen Intel Core i7-12700F, 1 CPU, 20 logical and 12 physical cores  
.NET SDK=8.0.100-preview.3.23178.7  
  [Host]     : .NET 8.0.0 (8.0.23.17408), X64 RyuJIT AVX2  
  DefaultJob : .NET 8.0.0 (8.0.23.17408), X64 RyuJIT AVX2  

Method	Mean	Error	StdDev	Code Size
Lerp1	18.78 ns	0.406 ns	0.965 ns	148 B
Lerp2	19.08 ns	0.448 ns	1.285 ns	168 B
Lerp3	18.91 ns	0.383 ns	0.974 ns	153 B
Lerp4	19.06 ns	0.411 ns	0.968 ns	173 B
Lerp5	19.29 ns	0.417 ns	0.843 ns	170 B

ほとんど変わりはないのですが、しいて言えば Lerp1 < Lerp2 = Lerp3 = Lerp4 < Lerp5 、という感じです。
文字通り誤差範囲内ですし、そのままだと測定不能になるレベルなのであまり気にしなくてもよさそうです。

擬似乱数生成との関係

ところで突然 Lerp の話を始めた理由はというと、実は擬似乱数生成と関係があるからです。
乱数生成において、 の範囲の乱数を返す NextDouble() をもとに、 の範囲に変換して NextDouble(min, max) とする処理はよくあります。
手書きで NextDouble() * (max - min) + min とした経験がある方も多いでしょう。

ここで役立つのが Lerp の式です。 の乱数を t として min, max 間を Lerp させることで、この範囲変換を行うことができます。

この範囲変換を行うにあたり、どの式にしたら一番いい感じになるのか（雑）、と思い立ったのが本記事のきっかけです。

おわりに

「全部同じじゃないですか～」「違いますよ～」といった幻聴がする気がしますが、実際どれも微妙に異なります。

精度を求めるなら fma(t, y - x, x) 、 t == 1 の場合や極端な値を与えたときの挙動を考えるなら fma(t, y, (1.0 - t) * x) が良い、でしょうか…？

おすすめの式はこれとか、実はこっちの式のほうがいいとか、有識者の方がいらっしゃいましたらコメント等いただけると助かります。

Unity で C ネイティブプラグインから __m128(i/d) 値を送受信する方法

はじめに

C プラグイン側にこういう関数があって、これを Unity (C#) 側から呼びたいとします。

__declspec(dllexport) __m128d nextDouble2(rng_t *rng, __m128d min, __m128d max)  
{  
    // do something...  
}  

なお、 __m128d は double 型を 2 つつなげた SIMD 用の型です。 float 4 つの __m128 と uint64_t 2 つの __m128i もあります。

Burst のヘルプ を見ると、

For all DllImport and internal calls, you can only use the following types as parameter or return types:
...
Unity.Burst.Intrinsics.v128

とあるので、 Burst を入れて v128 型で受ければよさそうに思えます。
v128 は 128 bit の union っぽい構造体で、例えば v128.ULong0 で下位 64 bit の ulong 値を取得できたりします。

というわけで C# 側で DllImport を書きましょう:

// not working :(  
[DllImport("DllName")]  
private static extern v128 nextDouble2(IntPtr rng, v128 min2, v128 max2);  

ところがこれは動きません。
呼び出すとランダムな値が返ったり、 Unity がフリーズしたりします。

解決策

前述のヘルプの下のほうを見ると、

Passing structs by value isn't supported; you need to pass them through a pointer or reference.

というわけで、サポートされていません。
いや私もこの部分は読んでいたのですが、まさか v128 も対象だとは思わず…
だって v128 (__m128) はビルトイン型みたいなものじゃないですか…

というわけで、呼べるようにするためにはポインタか ref 参照にしなければなりません。
ポインタにすると unsafe になって面倒なので、 ref でやります。

まずは C プラグイン側にラッパーを生やします。

__declspec(dllexport) void nextDouble2ForUnity  
    (rng_t* rng, __m128d* min2, __m128d* max2, __m128d* output) {  
    *output = nextDouble2(rng, *min2, *max2);  
}  

そして、 C# 側からこのようにして参照します。

// it works fine :)  
[DllImport("DllName")]  
private static extern void nextDouble2ForUnity  
    (IntPtr rng, in v128 min2, in v128 max2, out v128 output);  

ref の代わりに in と out を使うこともできるので、適宜使い分けます。

最後に C# 側のラッパーを書けば完成です。

private static v128 nextDouble2(IntPtr rng, v128 min2, v128 max2)  
{  
    nextDouble2ForUnity(rng, min2, max2, out var result);  
    return result;  
}  

おわりに

ググっても何も出てこないのでさっぱりわからず、 n 週間塩漬けにしていました…
私のようにうっかり勘違いをする人がいるかもしれないので残しておきます…

詳説 Ziggurat 法 ～ 正規分布・指数分布乱数の高速生成

はじめに

正規分布 や 指数分布 に従う擬似乱数を生成する際、 Ziggurat 法 を用いると高速に生成することができます。
本稿では、この Ziggurat 法の仕組みや実装について、できる限り詳しく解説していきたいと思います。

Ziggurat 法以外の手法について

早速脱線しますが、 Ziggurat 法以外にも正規分布・指数分布に従う乱数を生成する方法があります。
これらは実装が簡単であり、お世話になった方も多いかと思います。参考として挙げていきましょう。

なお、通常の System.Random では機能に不足があるので、以下のクラスを継承する擬似乱数生成器があるものとします。

public abstract class RandomBase  
{  
    /// <summary>  
    /// 64bit; [0, 2^64) の擬似乱数を生成   
    /// </summary>  
    public abstract ulong Next();  
}  

Box-Muller 法

Box-Muller 法 は、正規分布に従う乱数を 2 つ生成します。
C# では、以下のように実装します。

public static (double x, double y) BoxMuller<TRng>(TRng rng)  
    where TRng : notnull, RandomBase  
{  
    // [0, 1)  
    double u0 = (rng.Next() >> 11) * (1.0 / (1ul << 53));  
    double u1 = (rng.Next() >> 11) * (1.0 / (1ul << 53));  
  
    double r = Math.Sqrt(-2 * Math.Log(1 - u0));  
    (double sin, double cos) = Math.SinCos(u1 * (Math.PI * 2));  
  
    return (r * cos, r * sin);  
}  

念のため補足しておくと、 (rng.Next() >> 11) * (1.0 / (1ul << 53)) は ulong 型の乱数 から double 型の乱数 に変換するコードです。 double 型の精度は 53 bit のためこのようになっています。

Box-Muller 法は乱数の消費個数が固定なので (2 回の Next() から必ず 2 個得られる) 、乱数のストリームにフィルタをかける感じで使うなど、特定のシチュエーションでは便利です。

また余談ですが、 Wikipedia に載っている SVG が分かりやすくて素晴らしいので、一度触ってみることをおすすめします。

Polar 法

Polar 法 も、正規分布に従う乱数を 2 つ生成します。

public static (double x, double y) PolarMethod<TRng>(TRng rng)  
    where TRng : notnull, RandomBase  
{  
    double x, y, s;  
    do  
    {  
        // [-1, 1)  
        x = ((long)rng.Next() >> 10) * (1.0 / (1ul << 53));  
        y = ((long)rng.Next() >> 10) * (1.0 / (1ul << 53));  
  
        s = x * x + y * y;  
    } while (s >= 1 || s == 0);  
  
    double r = Math.Sqrt(-2 * Math.Log(s) / s);  
  
    return (x * r, y * r);  
}  

さらに念のため、 ((long)rng.Next() >> 10) * (1.0 / (1ul << 53)) は double 型の に変換するコードです。算術右シフトを活用して符号を設定しています。

Polar 法は乱数の消費個数がランダムです。
しかし、Box-Muller 法と比べると SinCos が不要なため、 一般に Polar 法 のほうが速くなります。

逆関数法 (指数分布)

指数分布に従う乱数を得るには、逆関数法が一番実装が簡単だと思います。

public static double Inverse<T>(T rng)  
    where T : notnull, RandomBase  
{  
    // [0, 1)  
    double u0 = (rng.Next() >> 11) * (1.0 / (1ul << 53));  
  
    return -Math.Log(1 - u0);  
}  

たったこれだけで済みます。
この式は、指数分布の累積分布関数 の逆関数 からきています。

ところで、 u0 は なので単に -Math.Log(u0) としてもほとんど違いはないような気もしますが、こうすると u0 == 0 の場合に ∞ を返してしまい、問題の種になるので素直に 1 - u0 としています。

Ziggurat 法

以上に挙げた手法では、 Log ・ Sqrt ・ SinCos といった数学関数を使用していました。
これらは四則演算よりは重い処理になるため、できる限り使用を回避したいところです。

Ziggurat 法では、高確率でこれらの数学関数を使用せずに正規分布・指数分布に従う乱数を生成することができます。
まずは概念から説明します。
なお、一旦正規分布を中心に説明を進めていきます。

概念

Ziggurat 法の原著論文は こちら を参照してください。 *1

Ziggurat 法は棄却採択法の一種です。
棄却採択法とは、欲しい分布自体の計算が難しい・時間がかかる場合に、より大きくて単純に計算できる分布・範囲に対してランダムに点を打ち、それが欲しい分布の範囲内なら採択・範囲外なら棄却してもう一回やる、という手順から乱数を得る手法です。

Ziggurat 法は、欲しい分布の上に面積が等しい長方形を敷き詰め、その長方形上にランダムに点を打ち、それが欲しい分布の中であれば採用する、という手法です。

これは正規分布に対して 8 レイヤーの Ziggurat を置いた例です。 8 つの赤い四角形 (最下段のみ右端へ続く "尾" を含みます) はすべて同じ面積ですので、生成が簡単な一様分布に従う乱数で選択することができます。
四角形を選んだら、さらにその内部にランダムに点を打ち、それが分布より下に (中に) なればその x 座標を採択 (return して) 、そうでなければ棄却 (やり直し)、という手順で生成します。

ざっくりした生成手順としては、

1. どのレイヤーにするか選ぶ
2. レイヤー内にランダムに点を打ち、長方形内部か、縁のエリアかの判定
   1. 長方形内部なら、その x 座標を返す
   2. 尾の部分なら、専用の計算をして返す
   3. 縁のエリアなら、確率密度関数の下にあるか判定
      1. 収まれば、その x 座標を返す
      2. 収まらなければ、最初からやり直す

という感じです。

ちなみに、 Ziggurat 法は対象の確率密度関数が で単調減少であるか、 軸に対称であれば適用可能です。
今回挙げた正規分布や指数分布以外にも、 コーシー分布 や ラプラス分布 などにも適用できます。

実装

テーブル生成

まずは各レイヤーの座標の定数テーブルを生成する必要があります。
これは分布に対して固定なので、プログラムの最初に生成するか、なんならすべて定数として埋め込んでもよいです。

/// <summary>  
/// Ziggurat 法のテーブル生成  
/// </summary>  
/// <param name="tableSize">128 or 256</param>  
/// <param name="probabilityDensityFunction">確率密度関数</param>  
/// <param name="cumulativeDistributionFunction">累積分布関数</param>  
/// <param name="inverseProbabilityDensityFunction">確率密度関数の逆関数</param>  
/// <returns>  
/// * 各レイヤー右上角の X 座標  
/// * 各レイヤー右上角の Y 座標  
/// * 各レイヤーの (確率密度関数と重ならない四角形の面積 / 全体の面積)  
/// </returns>  
public static (double[] rectangleX, double[] rectangleY, double[] insideRectangleRatio)  
    CreateZigguratTable(int tableSize,  
    Func<double, double> probabilityDensityFunction,  
    Func<double, double> cumulativeDistributionFunction,  
    Func<double, double> inverseProbabilityDensityFunction)  
{  
    double CalculateArea(double x0)  
    {  
        return x0 * probabilityDensityFunction(x0) + cumulativeDistributionFunction(double.PositiveInfinity) - cumulativeDistributionFunction(x0);  
    }  
  
    // 一番下のレイヤーの右上角の x 座標 ref. dn  
    double x0;  
    // 1 つのレイヤーの面積 ref. vn  
    double area;  
    {  
        // 誤差関数  
        double CalculateError(double x0)  
        {  
            double area = CalculateArea(x0);  
  
            double xi = x0;  
            for (int i = 1; i < tableSize; i++)  
            {  
                xi = inverseProbabilityDensityFunction(area / xi + probabilityDensityFunction(xi));  
  
                // 計算失敗時に負になる場合がある  
                if (xi < 0)  
                    return double.NegativeInfinity;  
            }  
  
            // 面積が足りない場合などに NaN になる場合がある  
            return double.IsNaN(xi) ? double.NegativeInfinity : xi;  
        }  
  
        // 探索を開始する最小値/最大値 (ヒューリスティック; 計算に失敗した場合は適宜変更する)  
        double left = 2;  
        double right = 10;  
  
        // 二分法で x0 を求める  
        while (true)  
        {  
            double mid = left * 0.5 + right * 0.5;  
  
            if (left == mid && mid < right)  
                mid = Math.BitIncrement(left);  
            else if (left < mid && mid == right)  
                mid = Math.BitDecrement(right);  
            if (left == mid || right == mid)  
                break;  
  
            double leftError = CalculateError(left);  
            double midError = CalculateError(mid);  
            double rightError = CalculateError(right);  
  
            if (Math.Sign(leftError) != Math.Sign(midError))  
                right = mid;  
            else if (Math.Sign(rightError) != Math.Sign(midError))  
                left = mid;  
            else  
                throw new InvalidOperationException("something wrong. should change left/right");  
        }  
  
        x0 = right;  
        area = CalculateArea(x0);  
    }  
  
  
    // 各レイヤー右上角の X 座標 ref. wn  
    var rectangleX = new double[tableSize];  
  
    // 各レイヤー右上角の Y 座標 ref. fn  
    var rectangleY = new double[tableSize];  
  
    // 各レイヤーの (確率密度関数と重ならない四角形の面積 / 全体の面積) ref. kn  
    var insideRectangleRatio = new double[tableSize];  
  
  
    // rectangleX の設定  
    // [0]: 一番下のレイヤーを長方形と見立てたときの横幅 (insideRectangleRatio[0] の設定のため)  
    // [1]: 一番下のレイヤー内の長方形の右端  
    rectangleX[0] = area / probabilityDensityFunction(x0);  
    rectangleX[1] = x0;  
    for (int i = 2; i < rectangleX.Length; i++)  
    {  
        double xi = inverseProbabilityDensityFunction(area / rectangleX[i - 1] + probabilityDensityFunction(rectangleX[i - 1]));  
        rectangleX[i] = xi;  
  
        if (xi < 0)  
            throw new InvalidOperationException("negative x; parameter may be wrong");  
    }  
  
    // rectangleY の設定  
    // [^1]: 頂点の y 座標  
    for (int i = 0; i < rectangleY.Length - 1; i++)  
    {  
        rectangleY[i] = probabilityDensityFunction(rectangleX[i + 1]);  
    }  
    rectangleY[^1] = probabilityDensityFunction(0);     // == 1;  
  
    // insideRectangleRatio の設定  
    // [^1]: 一番上のレイヤーはすべて範囲外なので 0  
    for (int i = 0; i < insideRectangleRatio.Length - 1; i++)  
    {  
        insideRectangleRatio[i] = rectangleX[i + 1] / rectangleX[i];  
    }  
    insideRectangleRatio[^1] = 0;  
  
  
    return (rectangleX, rectangleY, insideRectangleRatio);  
}  

与える引数はこんな感じです。

• tableSize : 定数テーブルのサイズ。通常 128 か 256
  • 2 冪にすることで高速に生成できるようにします。詳しくは後述
  • 一般に 256 のほうが速いですが、もちろん定数テーブルのサイズが大きくなります
• probabilityDensityFunction : 確率密度関数
  • 正規分布なら x => Math.Exp(-x * x / 2)
    • 正規化定数 (分布全体の面積を 1 にするために掛けられている定数) の は無視してかまいません
  • 指数分布なら x => Math.Exp(-x)
• cumulativeDistributionFunction : 累積分布関数
  • 正規分布なら Math.Sqrt(2 * Math.PI) * 0.5 * (1 + Erf(x / Math.Sqrt(2)))
    • Erf() は誤差関数 (後述)
  • 指数分布なら 1 - Math.Exp(-x)
• inverseProbabilityDensityFunction : 確率密度関数の逆関数
  • 正規分布なら y => Math.Sqrt(-2 * Math.Log(y))
  • 指数分布なら y => -Math.Log(y)

Erf() は 誤差関数 です。
C# の Math にはないので、自前でテイラー展開して実装するか Math.NET Numerics などから使うかしてください。

public static double Erf(double real)  
{  
    if (double.IsInfinity(real))  
        return Math.Sign(real);         // ±∞ => ±1  
  
    double erf = 2 / Math.Sqrt(Math.PI);  
    double sum = 0;  
    double factorial = 1;  
    for (int i = 0; true; i++)  
    {  
        double element = (((i & 1) != 0 ? -1 : 1) * Math.Pow(real, 2 * i + 1)) /  
            (factorial * (2 * i + 1));  
  
        if (sum + element == sum)  
            break;  
  
        sum += element;  
        factorial *= i + 1;  
    }  
    return sum * erf;  
}  

area と x0 の計算手法

area (各レイヤーの面積) や x0 (一番下のレイヤーの右角 x 座標) といった定数の計算手法について説明します。

tableSize が 256 のとき、面積 は以下の関係式を満たします:

は確率密度関数です。

2 番目の式を変形して、

とします。 は確率密度関数の逆関数で、 を満たすような関数です。
すると、 を決める → が決まる → が次々決まっていく → 最後の式

を満たしているかを確認する、といった手順で求めることはできそうです。

問題は最初の なのですが、どうやら数学的にうまいこと解くのは難しいらしく、求根アルゴリズム (二分法や割線法など) を使って求めることになります。

幸いにして はそこまで広い範囲にはないことが予想されるので (せいぜい 程度でしょう) 、今回は 二分法 を使って解きます。

二分法は、以下のような手順で となる を計算するアルゴリズムです。

1. 探索の最小値 left と最大値 right を決める
2. 中点 mid を求める
3. 符号に合わせて left か right を更新する
   1. sign(f(left)) != sign(f(mid)) なら right = mid
   2. sign(f(right)) != sign(f(mid)) なら left = mid
4. 以上を left と right が十分に狭い範囲になるまで繰り返す
5. left か right を目的の値として返す

area と x0 は、 tableSize に合わせて次のようになるはずです。

分布	tableSize	area	x0
正規分布	128	0.0099125630353356087	3.4426198558966847
正規分布	256	0.0049286732339721695	3.6541528853613281
指数分布	128	0.0079732295395533725	6.8983151166156444
指数分布	256	0.0039496598225815527	7.6971174701310288

定数テーブルの中身

CreateZigguratTable() で何が計算できるかを図で見てみましょう。
尾のあたりを拡大して見てみます:

rectangleX と rectangleY には各レイヤーの角の座標、 insideRectangleRatio にはオレンジ色のエリア (確率密度関数と交わる部分を含まない面積の割合) の情報が入っています。
rectangleX[0] は図にありませんが、これには別の値、具体的にはレイヤー 0 (一番下)を長方形と見立てたときの横幅が入っています。なぜこうするのかというと、 insideRectangleRatio[0] で (レイヤー 0 の長方形部分の面積 / レイヤー 0 全体の面積) を計算するためです。

頭のあたりを拡大するとこんな感じです。
(念のため、インデックスの [^i] は最後から i 番目 (最後の要素は [^1]) を表します。)

一番上のレイヤーの insideRectangleRatio[^1] は 0 になっていることに注意してください。

分布に従う乱数の生成

さて、テーブルの生成が終わったら、本題の乱数生成に入ります。

// CreateZigguratTable で生成した定数テーブル  
private static readonly double[] ZigguratNormalX = new double[] { ... };  
private static readonly double[] ZigguratNormalY = new double[] { ... };  
private static readonly double[] ZigguratNormalInsideRectangleRatio = new double[] { ... };  
  
public static double ZigguratNormal<TRng>(TRng rng)  
    where TRng : notnull, RandomBase  
{  
    const int TableSize = 256;  
  
    while (true)  
    {  
        long u1 = (long)rng.Next();  
        // [-1, 1)  
        double u1Double = (u1 >> 10) * (1.0 / (1ul << 53));  
        int rectangleIndex = (int)(u1 & (TableSize - 1));  
  
        // 確率密度関数の範囲内にある四角形なら、そのまま返す  
        if (Math.Abs(u1Double) < ZigguratNormalInsideRectangleRatio[rectangleIndex])  
            return u1Double * ZigguratNormalX[rectangleIndex];  
  
        if (rectangleIndex == 0)  
        {  
            // 一番下のレイヤーの尾の部分  
            double s, t;  
            double x0 = ZigguratNormalX[1];  
            do  
            {  
                s = -Math.Log(1 - (rng.Next() >> 11) * (1.0 / (1ul << 53))) / x0;  
                t = -Math.Log(1 - (rng.Next() >> 11) * (1.0 / (1ul << 53)));  
            } while (s * s > t + t);  
  
            return (x0 + s) * (u1 < 0 ? -1 : 1);  
        }  
        else  
        {  
            // それ以外の端の部分  
            double x = u1Double * ZigguratNormalX[rectangleIndex];  
            double y = (rng.Next() >> 11) * (1.0 / (1ul << 53));  
  
            if (ZigguratNormalY[rectangleIndex - 1] + y * (ZigguratNormalY[rectangleIndex] - ZigguratNormalY[rectangleIndex - 1]) <= Math.Exp(-0.5 * x * x))  
                return x;  
        }  
    }  
}  

順を追って説明していきましょう。

テーブル選択・四角形の中にあるか判定

long u1 = (long)rng.Next();  
// [-1, 1)  
double u1Double = (u1 >> 10) * (1.0 / (1ul << 53));  
int rectangleIndex = (int)(u1 & (TableSize - 1));  

まずは一様分布に従う乱数を生成して、 u1Double を -1.0 <= u1Double < 1.0 の範囲で、どのレイヤーを選ぶかを決める rectangleIndex を 0 <= rectangleIndex < tableSize の範囲で生成します。
テーブル生成にて tableSize を 2 冪にしておいたのは、ここで rectangleIndex をビットマスクだけで生成できるようにするためです。

// 確率密度関数の範囲内にある四角形なら、そのまま返す  
if (Math.Abs(u1Double) < ZigguratNormalInsideRectangleRatio[rectangleIndex])  
    return u1Double * ZigguratNormalX[rectangleIndex];  

次に、選択されたレイヤーの insideRectangleRatio[rectangleIndex] (オレンジの領域) 内に u1Double が入っているかを判定します。もし入っていればそれは確実に確率密度関数 (青線) の下にあるので、その x 座標、つまり u1Double * rectangleX[rectangleIndex] を返します。

ちなみに、正規分布で tableSize が 256 の場合は約 98.5 % の確率でここで採択されます。

そうでない場合は後続の処理に続きます。

尾の範囲外処理

// 一番下のレイヤーの尾の部分  
double s, t;  
double x0 = ZigguratNormalX[1];  
do  
{  
    // note: (rng.Next() >> 11) * (1.0 / (1ul << 53))  
    // は [0, 1) の乱数を生成する  
    s = -Math.Log(1 - (rng.Next() >> 11) * (1.0 / (1ul << 53))) / x0;  
    t = -Math.Log(1 - (rng.Next() >> 11) * (1.0 / (1ul << 53)));  
} while (s * s > t + t);  
  
return (x0 + s) * (u1 < 0 ? -1 : 1);  

選ばれたのが一番下のレイヤーで、かつ範囲外だった場合は一工夫が必要です。
正規分布では、上記のコードで尾の部分の乱数を生成することができます。

……どうして?

説明していきましょう。こちらについては、文献 *2 を参考にしつつ解説します。

まず、正規分布の確率密度関数から正規化定数を除いた と、一番下のレイヤーの右角 (rectangleX[1]) である について、

となるので、 は常に正規分布の確率密度関数より大きくなります。したがって棄却採択法の関数として使えます。
(繰り返しになりますが、計算が簡単でより大きな分布を持つ関数で覆ってしまって、その内部に点を打ち、それが計算が複雑で小さい分布 (ここでは正規分布) の中にあれば OK 、というのが棄却採択法です。)

確率密度を求めるため、 の範囲で積分すると、

になりますので、これで当初の式を割って正規化すると

を得ます。これが確率密度関数になります。
この累積分布関数は、

となります。この逆関数を求めましょう。

注意すべきポイントとしては、 とすると なので となってしまうので、そうならないように の範囲を として にしています。

ここで、 の範囲の一様分布乱数 を用意して、 に を掛けた がもともとの を下回れば分布の中に入るので採択、そうでなければ棄却、としましょう。

先ほど求めた逆関数 より、 ですから、

を得ます。
さらに とおけば、

ここまで来れば見覚えのある式になったのではないでしょうか?
コードをもう一度見直してみましょう。

// 一番下のレイヤーの尾の部分  
double s, t;  
double x0 = ZigguratNormalX[1];  
do  
{  
    s = -Math.Log(1 - (rng.Next() >> 11) * (1.0 / (1ul << 53))) / x0;  
    t = -Math.Log(1 - (rng.Next() >> 11) * (1.0 / (1ul << 53)));  
} while (s * s > t + t); // 棄却されたらやり直しなので、不等号が逆  
  
return (x0 + s) * (u1 < 0 ? -1 : 1);  

このような処理になった理由が分かったかと思います。

また、 こちらのブログ にも別のアプローチも含めて詳しく述べられています。

---

また、指数分布ではこのように処理します。

// 一番下のレイヤー (尾の部分)  
double x0 = ZigguratExponentialX[1];  
return x0 + ZigguratExponential(rng);  

指数分布の尾の部分の処理は、単に を足して指数分布と同じ手法を再帰的に適用するだけでよいです。
なぜかを説明します。

より右側の部分の面積は、

です。次に、 以上の部分の累積分布関数を求めます。
通常時の累積分布関数の始点を にしたうえで面積で割ることで、

となりますので、その逆関数は、

となります。これは、通常時の累積分布関数の逆関数 に を足したものになります。

というわけで、 を足して指数分布に従う乱数を再帰的に生成すればよいことになります。
(再帰すると遅いので、実際はループで処理します。)

それ以外の範囲外処理

さて、尾以外の部分で範囲外となった場合は、以下のように処理していました。

// それ以外の端の部分  
double x = u1Double * ZigguratNormalX[rectangleIndex];  
double y = (rng.Next() >> 11) * (1.0 / (1ul << 53));  
  
if (ZigguratNormalY[rectangleIndex - 1] + y * (ZigguratNormalY[rectangleIndex] - ZigguratNormalY[rectangleIndex - 1]) <= Math.Exp(-0.5 * x * x))  
    return x;  

位置関係を図にするとこんな感じになります。

こういう感じで座標を計算して、それが確率密度関数よりも低ければ採択、そうでなければ棄却してやり直します。

高速化ポイント

オリジナルの実装では、 insideRectangleRatio を整数型にして直接 u1 (整数乱数) と比較することで、浮動小数点数型同士ではなく整数型同士の比較にして高速化を図っているようです。

また、 rectangleX に などをあらかじめ掛けておくことで、 u1Double のシフト・定数との掛け算を省いて直接 u1 * rectangleX[rectangleIndex] として返す、といった手法もあります。

今回は分かりやすさを重視したため以上の手法は使っていませんが、高速化の余地はあります。

改良 Ziggurat 法

Ziggurat 法には改良されたバージョンがあります。 *3
改良 Ziggurat 法は Java の nextGaussian() ・ nextExponential() の実装でも 使われています。

こちらの紹介もしていきたいと思います。
なんなら今までが前置きでここから本題です。

概念 ～ 改良されたポイント

Ziggurat の配置について

改良 Ziggurat 法の Ziggurat を見てみましょう。

これは、正規分布に対して 8 レイヤーの Ziggurat を置いた例です。
……はい。 8 レイヤーです。でも 5 レイヤーしかありませんし、確率密度関数からはみ出していません。
これでは破綻しているのでは? と思われるかもしれませんので、説明を続けます。

残りの 3 レイヤーはどこにあるのかというと、四角形の外と確率密度関数に挟まれた端の (オレンジ色の) 部分です。
実は端の部分の面積をすべて足すと、ちょうどレイヤー 3 枚分になります。
これによって、「同じ面積のレイヤーをランダムに (一様分布の乱数で) 選択する」機能を維持することができます。
また、それ以外の四角形のレイヤーが選択された場合に、端の処理をせずに即採択することができます。厄介だった端の処理を 3 枚分に押し込めたイメージです。
正規分布では、全体の約 98.8 % が四角形のレイヤーで採択されます。よく通るパスで端の処理を省略できるため、高速化につながります。

それでは実際に端の 3 枚分が選択された場合はどのようにして生成するのかというと、

1. 各レイヤーの端の部分の面積を重みにした重み付き抽選を行い、レイヤーインデックス j を得る
2. Ziggurat 法と同様に、レイヤー j における端の部分の処理を行う
   1. レイヤー 0 (一番下) の場合は、尾から生成する処理を行う
   2. それ以外のレイヤーの場合は、確率密度関数より下になるかを判定して棄却/採択する

という感じで、レイヤーを再抽選したあとは Ziggurat 法で端を引いた時と似たような処理になります。

最初に重み付き抽選を行うのは、もちろん各レイヤーの端の部分の面積がそれぞれ異なるからです。
先ほどの図を見ていただければわかる通り、一番上のレイヤーの端は巨大で、その下のレイヤーの端と比べると一目瞭然でしょう。
そのため、その面積に従って重み付き抽選をする必要があります。

重み付き抽選を安直にやってしまうと時間がかかるので、 Walker's Alias method を利用します。

Walker's Alias method

Walker's Alias method は、事前にテーブル生成をすることで重み付き抽選を高速に行うアルゴリズムです。 *4

通常、重み付き抽選には 、重み順にソートして二分探索するなどしても はかかりますが、 Alias method を利用すると抽選は で行うことができます。驚きの速さです。

その分事前に かかるテーブル生成をする必要があるため、重みを場合によって変えたり、重みテーブルを 1 回しか使わなかったりする場合では安直な方法より遅くなりがちです。重みを変えずに何度も抽選する場合に有効なアルゴリズムであるといえます。
今回適用したい確率分布は絶対に重みが変わらないため事前にテーブル生成を済ませておくことができ、まさにおあつらえ向きな手法です。

実例を挙げると、前述した各レイヤーの端部分の面積はこのようになっています。
一番上のレイヤーが重い (大きい) ことがよくわかります。

これを、

• 全てのレイヤーが同じ重さになるように、 8 つのレイヤーに切り分ける
• 1 つのレイヤーに 2 種類まで入れる
• 2 種類のうち 1 つは元々あったもの

というルールでいい感じに分割すると、以下のようになります。

端に振ってある数字が元々のレイヤーのインデックス (一番下が 0 番目) です。
大きかった 5 番目 (一番上) がまんべんなく振られていますね。
ちなみに、上 3 つの左側には重さ 0 の青いほうがあるものとします。

この振り分けが終わると、以下の手順で重み付き抽選をすることができます。

1. 乱数を生成し、レイヤーをランダムに選ぶ。
2. 0 ～ 1 の乱数を生成し、閾値未満なら元々のインデックス (左の青いほう) 、閾値以上なら 2 番目のインデックス (右の赤いほう) を返す。

これはテーブルを毎回走査する必要がなく、乱数消費も 2 個固定です。したがって で重み付き抽選が実現されます。

このテーブルが持つべき情報は、「左側の重み」「右側のインデックス」です。
右側の重みは (1 - 左側の重み) になりますし、左側のインデックスはそのレイヤーのインデックスと同じにしてあるためです。

なお、詳しいテーブル作成手法などについては、 Wikipedia や こちらの記事 が参考になります。

端の部分の最適化

さらに最適化を進めます。

端の部分の判定を行うとき、今までは単に「確率密度関数より下に行くか?」を見ていましたが、ここにも最適化の余地があります。

端の部分は大体三角形をしています。中でも、「ふくらんだ三角形」「へこんだ三角形」「ふくらみとへこみのある三角形」に分類することができます。

ふくらんだ三角形

一番上のレイヤーを例に挙げて説明しましょう。

まず、この四角形の領域に斜線を引いて三角形を作りましょう。
すると、ふくらんだ三角形の場合、左下のオレンジの三角形の部分は常に青線の確率密度関数より下になることがわかります。
したがって、ランダムに打った点がオレンジの領域に入った場合は、確率密度関数の計算を省略して採択することができます。
確率密度関数の計算では Math.Exp(-0.5 * x * x) や Math.Exp(-x) など相対的に重い指数関数を使いますので、これを省略できるとスピードアップにつながります。

また、一番ふくらんでいるところに接するように、同じ傾きの線を引きます。
すると、右上の水色の三角形の部分は常に青線の確率密度関数より上になり、同様に確率密度関数の計算を省略して棄却することができます。

なお、オレンジと水色の三角形に挟まれた部分はどうしようもないので、今までと同様に確率密度関数の計算を行ってそれより下かを判定します。

へこんだ三角形

下から二番目のレイヤーを例に挙げて説明しましょう。

上と同様に、一番へこんでいるところに接するように同じ傾きの線を引くと、左下のオレンジ色の部分は必ず採択される領域になります。

また、左上の部分はそのままだと必ず棄却される領域になりますが、もったいないので 180°回転させる (または上下左右を反転させる) ことで再抽選を行います。
再抽選でも同様に、左下の緑色の部分は必ず採択、残った水色の部分は諦めて確率密度関数との計算を行うようにします。

なお、指数分布ではふくらみはないので、全ての領域がへこんだ三角形になります。

ふくらみとへこみのある三角形

ふくらみとへこみの中間地点がこれにあたります。
ここでも同様に、

• 一番へこんでいるところに接するように同じ傾きの線を引き、それより下は採択
• 一番ふくらんでいるところに接するように同じ傾きの線を引き、それより上は棄却
• その中間では通常通り、確率密度関数との判定

とします。

尾の生成の最適化

さて、一番下のレイヤーの尾の部分についても考えてみましょう。

以前は尾の部分の生成ではこういう感じのコードになっていました。

// 一番下のレイヤーの尾の部分  
double s, t;  
double x0 = ZigguratNormalX[1];  
do  
{  
    s = -Math.Log(1 - (rng.Next() >> 11) * (1.0 / (1ul << 53))) / x0;  
    t = -Math.Log(1 - (rng.Next() >> 11) * (1.0 / (1ul << 53)));  
} while (s * s > t + t);  
  
return (x0 + s) * (u1 < 0 ? -1 : 1);  

ここで、 -Math.Log(...) の部分に着目してみましょう。
これは、冒頭で述べた「逆関数法による指数分布に従う乱数生成」と同じになります。
したがって、 (改良) Ziggurat 法による生成のほうが速いのであれば、それに置き換えることができます。

// 最下層レイヤー: 尾から生成する  
double x0 = ModifiedZigguratNormal128X[0] * (1ul << 63);  
double s, t;  
do  
{  
    s = ModifiedZigguratExponential(rng) / x0;  
    t = ModifiedZigguratExponential(rng);  
} while (s * s > t + t);  
return (x0 + s) * (u1 < 0 ? -1 : 1);  

実装

以上を踏まえて、実装に移っていきましょう。
まずは定数テーブル生成からやります。

定数テーブル生成

改良 Ziggurat 法においては、面積 () は単純に (全域の面積 / tableSize) で求めることができます。はみ出す部分がないのでここの計算は楽です。

今回の正規分布での関係式は、

となります。なお、 は累積分布関数です。
まずは最後の式から を求めることができます。
そして 2 番目の式から、

となるので、 から を求めることができます。
したがって、 (上から 2 番目のレイヤーの四角形右端) を決める → が決まっていく → が決まる → 最初の式 を満たすかを確認する、といった手順で進められます。

Ziggurat 法とは決定の順序が異なり、上から順 (インデックスが大きな順) なことに注意してください。
また、指数分布の場合は端の面積が 3 レイヤーだと足りず 4 レイヤー分になるので、 3 (252) を 4 (251) に読み替えてください。

/// <summary>  
/// 改良 Ziggurat 法のテーブル生成  
/// </summary>  
/// <param name="tableSize">128 or 256</param>  
/// <param name="iMax">正規分布は <paramref name="tableSize"/>-3, 指数分布は <paramref name="tableSize"/>-4 </param>  
/// <param name="isSymmetric">正規分布は true, 指数分布は false</param>  
/// <param name="probabilityDensityFunction">確率密度関数</param>  
/// <param name="cumulativeDistributionFunction">累積分布関数</param>  
/// <param name="inversePdf">確率密度関数の逆関数</param>  
/// <param name="derivativePdf">確率密度関数の傾き (微分)</param>  
/// <returns>  
/// * X に係数 2^-n を掛けたもの; n=63(正規分布) or 64(指数分布)  
/// * Y に係数 2^-n を掛けたもの  
/// * 上位 56 bit が 1 つめの重み, 下位 8 bit が 2 つめのインデックスである Walker's Alias method のテーブル  
/// * 変曲点 (確率密度関数のへこみとふくらみの境界) のレイヤーのインデックス  
/// * 最大のふくらみの比  
/// * 最大のへこみの比  
/// </returns>  
public static (double[] x, double[] y, ulong[] aliasBox, int inflectionIndex, ulong maxConvexRate, ulong maxConcaveRate)  
    CreateModifiedZigguratTable(int tableSize, int iMax, bool isSymmetric,  
    Func<double, double> probabilityDensityFunction,  
    Func<double, double> cumulativeDistributionFunction,  
    Func<double, double> inversePdf,  
    Func<double, double> derivativePdf)  
{  
    // 1 レイヤーの面積  
    var area = (cumulativeDistributionFunction(double.PositiveInfinity) - cumulativeDistributionFunction(0)) / tableSize;  
  
    // 二分法  
    static double Bisection(Func<double, double> errorFunction, double initialLeft, double initialRight)  
    {  
        double left = initialLeft;  
        double right = initialRight;  
  
        while (true)  
        {  
            double mid = left * 0.5 + right * 0.5;  
  
            if (left == mid && mid < right)  
                mid = Math.BitIncrement(left);  
            else if (left < mid && mid == right)  
                mid = Math.BitDecrement(right);  
            if (left == mid || right == mid)  
                break;  
  
            double leftError = errorFunction(left);  
            double midError = errorFunction(mid);  
            double rightError = errorFunction(right);  
  
            if (Math.Sign(leftError) != Math.Sign(midError))  
                right = mid;  
            else if (Math.Sign(rightError) != Math.Sign(midError))  
                left = mid;  
            else  
                throw new InvalidOperationException("something wrong. should change left/right");  
        }  
  
        return left;  
    }  
  
  
    // xtop (x[iMax]) を二分法で求める  
    double xtop;  
    {  
        // 探索の最小/最大値。(ヒューリスティック; 動かないときは適宜変更する)  
        double xtopLeft = 0.01;  
        double xtopRight = 0.5 * (isSymmetric ? 1.0 : 0.4);  
  
        // 探索の誤差関数。最後のレイヤーの x 座標まで計算可能なら負の、不可能なら正の値を返す  
        double CalculateError(double xtop)  
        {  
            for (int i = iMax - 1; i >= 0; i--)  
            {  
                xtop = inversePdf(probabilityDensityFunction(xtop) - area / xtop);  
  
                if (double.IsNaN(xtop))  
                    return i + 1;  
            }  
            return -xtop;  
        }  
  
        xtop = Bisection(CalculateError, xtopLeft, xtopRight);  
    }  
  
  
    // x, y の設定  
    var x = new double[iMax + 1];  
    var y = new double[iMax + 1];  
  
    x[^1] = 0;  
    y[^1] = probabilityDensityFunction(0);  
    x[^2] = xtop;  
    y[^2] = probabilityDensityFunction(xtop);  
    for (int i = x.Length - 3; i >= 0; i--)  
    {  
        x[i] = inversePdf(probabilityDensityFunction(x[i + 1]) - area / x[i + 1]);  
        y[i] = probabilityDensityFunction(x[i]);  
  
        // verification  
        double actual = x[i + 1] * (y[i + 1] - y[i]);  
        if (Math.Abs(actual - area) > 1.0 / (1ul << 52) || double.IsNaN(actual))  
            throw new InvalidOperationException("area mismatch");  
    }  
  
  
    // 各レイヤー端の面積の計算  
    var edgeArea = new double[iMax + 1];  
    edgeArea[0] = cumulativeDistributionFunction(double.PositiveInfinity) - cumulativeDistributionFunction(x[0]);  
    for (int i = 1; i < edgeArea.Length; i++)  
    {  
        edgeArea[i] = (cumulativeDistributionFunction(x[i - 1]) - cumulativeDistributionFunction(x[i]))  
            - y[i - 1] * (x[i - 1] - x[i]);  
    }  
    // verification  
    {  
        double actual = edgeArea.Sum();  
        if (Math.Abs(actual - area * (tableSize - iMax)) > 1.0 / (1ul << 50))  
            throw new InvalidOperationException("edgeArea mismatch");  
    }  
  
  
    // Walker's Alias method によるテーブル作成  
    var firstWeights = new double[tableSize];  
    var secondIndexes = Enumerable.Range(0, tableSize).ToArray();  
    double edgeAverage = (tableSize - iMax) * area / tableSize;  
    {  
  
        for (int i = 0; i < edgeArea.Length; i++)  
            firstWeights[i] = edgeArea[i];  
  
        // 小さい順に並べ替える  
        var sorted = Enumerable.Range(0, tableSize).ToArray();  
        Array.Sort(sorted, Comparer<int>.Create((a, b) => firstWeights[a].CompareTo(firstWeights[b])));  
  
        // 最大の要素から最小の要素に分配していく  
        int small = 0;  
        int large = firstWeights.Length - 1;  
        while (small < large)  
        {  
            secondIndexes[sorted[small]] = sorted[large];  
  
            firstWeights[sorted[large]] -= edgeAverage - firstWeights[sorted[small]];  
  
            switch (firstWeights[sorted[large]].CompareTo(edgeAverage))  
            {  
                case > 0:  
                    small++;  
                    break;  
                case 0:  
                    small++;  
                    firstWeights[sorted[large]] = edgeAverage;  
                    secondIndexes[sorted[large]] = sorted[large];  
                    large--;  
                    break;  
                case < 0:  
                    (sorted[small], sorted[large]) = (sorted[large], sorted[small]);  
                    large--;  
                    break;  
            }  
        }  
  
    }  
    // firstWeights 正規化と検証  
    for (int i = 0; i < firstWeights.Length; i++)  
    {  
        if (firstWeights[i] - edgeAverage > 1.0 / (1ul << 52))  
            throw new InvalidOperationException("weight mismatch");  
  
        firstWeights[i] /= edgeAverage;  
    }  
  
  
    // 変曲点・ふくらみ/へこみの最大値の探索  
    int inflectionIndex = tableSize;  
    double maxConvexRatio = 0;  
    double maxConcaveRatio = 0;  
    for (int i = 1; i < x.Length; i++)  
    {  
        // > 0 ならふくらみ(convex), < 0 ならへこみ(concave)  
        double GradientDifference(double t)  
        {  
            double pdfGradient = derivativePdf(t);  
            double lineGradient = (y[i] - y[i - 1]) / (x[i] - x[i - 1]);  
            return pdfGradient - lineGradient;  
        }  
  
        double maxGradientPoint = Bisection(GradientDifference, x[i], x[i - 1]);  
  
  
        if (inflectionIndex == tableSize && GradientDifference(maxGradientPoint) > 0)  
        {  
            inflectionIndex = i;  
        }  
        {  
            double lineGradient = (y[i] - y[i - 1]) / (x[i] - x[i - 1]);  
            double diff = probabilityDensityFunction(maxGradientPoint) -  
                (lineGradient * maxGradientPoint + y[i] - lineGradient * x[i]);  
            double ratio = diff / (y[i] - y[i - 1]);  
  
            if (i < inflectionIndex)  
                maxConcaveRatio = Math.Max(-ratio, maxConcaveRatio);  
            else  
                maxConvexRatio = Math.Max(ratio, maxConvexRatio);  
        }  
    }  
  
  
    // x, y のスケール  
    for (int i = 0; i < x.Length; i++)  
    {  
        x[i] *= (isSymmetric ? 1.0 : 0.5) / (1ul << 63);  
        y[i] *= (isSymmetric ? 1.0 : 0.5) / (1ul << 63);  
    }  
  
    // Walker's Alias method のテーブルを一つにまとめる  
    var aliasBox = new ulong[tableSize];  
    for (int i = 0; i < aliasBox.Length; i++)  
    {  
        ulong ulongWeight = (ulong)(firstWeights[i] * (2.0 * (1ul << 63)));  
  
        // verification  
        if ((ulongWeight & ((ulong)tableSize - 1)) != 0)  
            throw new InvalidOperationException("weight may lose precision");  
  
        aliasBox[i] = ulongWeight ^ (ulong)secondIndexes[i];  
    }  
  
    return (x, y, aliasBox, inflectionIndex, (ulong)(maxConvexRatio * (2.0 * (1ul << 63))), (ulong)(maxConcaveRatio * (2.0 * (1ul << 63))));  
}  

以下、要所の説明をしていきます。

ふくらみ・へこみの計算

ふくらみ・へこみの計算では、「確率密度関数の傾き」が「四角形の左上と右下を結ぶ直線の傾き」が同じになる (または近づく) ポイントを探します。そこがふくらみ/へこみ (確率密度関数と直線との高さの差) が最大化されるポイントになります。
この図をもう一度見ていただくとわかりやすいかと思います。

ここで、ふくらみ・へこみは、四角形の高さとの割合として記録しておき、かつ全レイヤー中の最大値だけを取っておくようにします。これによって定数テーブルを増やさずに済みます。
原著論文曰く、これをテーブル化しても速度の向上につながらなかったそうです。 (Appendices D の 2 を参照)

X ・ Y のスケーリング

x と y に 1.0 / (1ul << 63) () を掛けていますが、これは後々 long 型の乱数 (範囲が になります) を直接掛けて の範囲にすることを想定しているためです。

double 型の仕様上、オーバーフロー・アンダーフローしない限り 2 冪の掛け算で精度が落ちることはありません。

Walker's Alias method のテーブルをまとめる

本来 Walker's Alias method のテーブルは「1 つめの要素の重み」「2 つめの要素のインデックス」の 2 つが必要ですが、これをまとめてしまおうという試みです。

重み firstWeights は正規化してあるので の範囲になっているのですが、これは double 型なのでだいたい 53 bit 程度の精度があります。これを ulong 型 (64 bit) に変換するのですが、 11 bit ほど余ることになります。
ここで、インデックス secondIndexes は なので 8 bit あれば足りますので、一緒の変数に入れてしまおう、という発想です。
メモリアクセスを減らすことができるので、多少は高速化が見込める……かもしれません。

分布に従う乱数の生成

ようやくテーブル生成が終わりました。さて、実際に乱数を生成するコードの説明に移っていきましょう。

// CreateModifiedZigguratTable で生成した定数テーブル  
private static readonly double[] ModifiedZigguratNormalX = { ... };  
private static readonly double[] ModifiedZigguratNormalY = { ... };  
private static readonly ulong[] ModifiedZigguratNormalAliasBox = { ... };  
  
[MethodImpl(MethodImplOptions.AggressiveInlining)]  
public static double ModifiedZigguratNormal<TRng>(TRng rng)  
    where TRng : notnull, RandomBase  
{  
    const int TableSize = 256;  
    // 四角形のレイヤー数  
    const int IMax = 253;  
  
    long u1 = (long)rng.Next();  
    int tableX = (int)(u1 & (TableSize - 1));  
  
    // 四角形のレイヤーなら x 座標を返す  
    if (tableX < IMax)  
    {  
        return u1 * ModifiedZigguratNormalX[tableX];  
    }  
  
    // 以降の処理は相対的に重いため別関数に分け、上の高速パスのインライン化を試みる  
    [MethodImpl(MethodImplOptions.NoInlining)]  
    static double Fallback(TRng rng, long u1)  
    {  
        // CreateModifiedZigguratTable で生成した定数  
        const int InflectionIndex = 204;  
        const ulong MaxConvexRate = 0x3efb83be6450cc00;  
        const ulong MaxConcaveRate = 0x151b6b6b7cd81f00;  
  
        [MethodImpl(MethodImplOptions.AggressiveInlining)]  
        static double SampleX(double x, int tableY) => Math.FusedMultiplyAdd(x, ModifiedZigguratNormalX[tableY - 1] - ModifiedZigguratNormalX[tableY], ModifiedZigguratNormalX[tableY] * (1ul << 63));  
        [MethodImpl(MethodImplOptions.AggressiveInlining)]  
        static double SampleY(double y, int tableY) => Math.FusedMultiplyAdd(y, ModifiedZigguratNormalY[tableY - 1] - ModifiedZigguratNormalY[tableY], ModifiedZigguratNormalY[tableY] * (1ul << 63));  
  
        // Walker's Alias method のテーブルを引く  
        ulong u2 = rng.Next();  
        int tableAlias = (int)(u2 & (TableSize - 1));  
        int tableY = (u2 >> 8) < (ModifiedZigguratNormalAliasBox[tableAlias] >> 8) ? tableAlias : (int)(ModifiedZigguratNormalAliasBox[tableAlias] & (TableSize - 1));  
  
        if (tableY == 0)  
        {  
            // 最下層レイヤー: 尾から生成する  
            double x0 = ModifiedZigguratNormalX[0] * (1ul << 63);  
            double s, t;  
            do  
            {  
                s = ModifiedZigguratExponential(rng) / x0;  
                t = ModifiedZigguratExponential(rng);  
            } while (s * s > t + t);  
            return (x0 + s) * (u1 < 0 ? -1 : 1);  
        }  
        else if (tableY < InflectionIndex)  
        {  
            // へこみレイヤー  
            ulong ux = (ulong)u1 << 1 >> 1;  
            ulong uy = rng.Next() >> 1;  
  
            double tx;  
  
            for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)  
            {  
                // 三角形の右上; 上下左右反転させる  
                if (uy < ux)  
                {  
                    ux = (1ul << 63) - 1 - ux;  
                    uy = (1ul << 63) - 1 - uy;  
                }  
  
                // 三角形の左下; 採択  
                if (uy >= ux + MaxConcaveRate)  
                {  
                    tx = SampleX(ux, tableY);  
                    break;  
                }  
  
                // 確率密度関数より下なら採択  
                tx = SampleX(ux, tableY);  
                double ty = SampleY(uy, tableY);  
  
                if (ty <= Math.Exp(-0.5 * tx * tx))  
                {  
                    break;  
                }  
            }  
            return tx * (u1 < 0 ? -1 : 1);  
        }  
        else if (tableY == InflectionIndex)  
        {  
            // へこみ＋ふくらみレイヤー  
            ulong ux = (ulong)u1 << 1 >> 1;  
            ulong uy = rng.Next() >> 1;  
  
            double tx;  
  
            for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)  
            {  
                // 三角形の左下; 採択  
                if (uy >= ux + MaxConcaveRate)  
                {  
                    tx = SampleX(ux, tableY);  
                    break;  
                }  
  
                // 三角形の右上; 棄却  
                if (uy + MaxConvexRate <= ux)  
                {  
                    continue;  
                }  
  
                // 確率密度関数より下なら採択  
                tx = SampleX(ux, tableY);  
                double ty = SampleY(uy, tableY);  
  
                if (ty <= Math.Exp(-0.5 * tx * tx))  
                {  
                    break;  
                }  
            }  
            return tx * (u1 < 0 ? -1 : 1);  
        }  
        else  
        {  
            // ふくらみレイヤー  
            ulong ux = (ulong)u1 << 1 >> 1;  
            ulong uy = rng.Next() >> 1;  
  
            double tx;  
  
            for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)  
            {  
                // 三角形の左下; 採択  
                if (uy >= ux)  
                {  
                    tx = SampleX(ux, tableY);  
                    break;  
                }  
  
                // 三角形の右上; 棄却  
                if (uy + MaxConvexRate <= ux)  
                {  
                    continue;  
                }  
  
                // 確率密度関数より下なら採択  
                tx = SampleX(ux, tableY);  
                double ty = SampleY(uy, tableY);  
  
                if (ty <= Math.Exp(-tx * tx / 2))  
                {  
                    break;  
                }  
            }  
            return tx * (u1 < 0 ? -1 : 1);  
        }  
    }  
    return Fallback(rng, u1);  
}  

大体の工夫は既に述べたとおりですが、一点あるとすればインライン化周りです。

インライン化されたメソッドは、メソッド内部のコードが呼び出し側に埋め込まれます。メソッド呼び出しのオーバーヘッドを省略できるため高速化につながります。 (++C++ を参照 。)
そのため MethodImplAttribute の AggressiveInlining で積極的にインライン化してほしい指定を行っています。

一方、 Fallback() ではわざわざそれを阻害する NoInline を指定しています。これは珍しいかもしれませんが、一応意図があります。
というのも、前述したように全体の約 98 % は Fallback() に到達する前に return します。
そして、 Fallback() には相対的に重く大きいコード (ループなど) が含まれています。したがって、 Fallback() のコードが展開された結果 ModifiedZigguratNormal() 本体のインライン化を阻害する結果にならないよう、あえてメソッド呼び出しとすること (NoInline) を指定しているというわけです。

パフォーマンス

ここまでに挙げた手法のパフォーマンスを比較してみましょう。
BenchmarkDotNet を使って計測しました。

ついでに、正規分布の改良 Ziggurat 法については「AliasBox を Weight と Index の 2 変数に分離した Separate」「ふくらみ・へこみの処理を省略してコードサイズを削った Simple」「tableSize = 128 の 128」についても試して、各手法の有効性を検証してみました。

BenchmarkDotNet=v0.13.5, OS=Windows 11 (10.0.22621.1555/22H2/2022Update/SunValley2)  
12th Gen Intel Core i7-12700F, 1 CPU, 20 logical and 12 physical cores  
.NET SDK=7.0.202  
  [Host]     : .NET 7.0.4 (7.0.423.11508), X64 RyuJIT AVX2  
  DefaultJob : .NET 7.0.4 (7.0.423.11508), X64 RyuJIT AVX2  

Method	Mean	Error	StdDev	Ratio
Inverse	6.7082 ns	0.1553 ns	0.3690 ns	1.00
ZigguratExponential	4.2090 ns	0.1120 ns	0.2157 ns	0.62
ModifiedZigguratExponential	5.5464 ns	0.1279 ns	0.1570 ns	0.79
BoxMuller	15.3776 ns	0.2602 ns	0.2307 ns	1.00
PolarMethod	13.0116 ns	0.2855 ns	0.3288 ns	0.85
ZigguratNormal	3.3800 ns	0.0815 ns	0.0762 ns	0.22
ModifiedZigguratNormal	0.9127 ns	0.0467 ns	0.0923 ns	0.06
ModifiedZigguratNormalSeparate	1.0547 ns	0.0484 ns	0.0967 ns	0.07
ModifiedZigguratNormalSimple	0.9783 ns	0.0473 ns	0.0750 ns	0.06
ModifiedZigguratNormal128	1.4022 ns	0.0557 ns	0.0704 ns	0.09

なお、 BoxMuller と Polar は一度に 2 個生成するので、実効の時間はこの半分として考えてください。

指数分布では、元々逆関数法もシンプルなためか速度の差は大きくありませんが、それでも改善しています。
また、なぜか改良前の Ziggurat 法のほうが速くなっています。原因については…コードサイズでしょうか…?

正規分布では、 Ziggurat 法が Box-Muller 法・ Polar 法に大きな差をつけており、改良 Ziggurat 法に至っては 1 ns を切っています。
また、 AliasBox で 1 変数にまとめたこと・ふくらみとへこみの処理をきちんと行うことによって、わずかではあるものの加速できていることも分かりました。
また、 tableSize を 128 にした場合 1.5 倍ほどかかるようになってしまいましたが、それでも Polar 法の 5 倍近く高速です。定数テーブルの削減を図りたい場合は採用の余地がありそうです。

おわりに

とても長くなってしまいましたが、 Ziggurat 法と改良 Ziggurat 法について解説しました。
あまり日本語の文献がないので難儀しましたが、この記事をまとめていくうえで私自身の理解も深めることができました。参考になれば幸いです。
もしわかりにくい点・怪しい点などあればご指摘ください。

Ziggurat 法は実装が相対的に複雑なため、ぱっと書けるようなものではないですが、定数テーブルさえ埋め込んでしまえばあとは比較的楽なので、活用できると良さそうです。

最後に、以上に示したコード群をまとめて置いておきます。ご活用ください。
本文では記載しなかった指数分布のコードも含んでいます。

Ziggurat / Modified Ziggurat method

*1:Marsaglia, G., & Tsang, W. W. (2000). The ziggurat method for generating random variables. Journal of statistical software, 5, 1-7.

*2:四辻哲章. (2010). 計算機シミュレーションのための確率分布乱数生成法. プレアデス出版.

*3:McFarland, C. D. (2016). A modified ziggurat algorithm for generating exponentially and normally distributed pseudorandom numbers. Journal of statistical computation and simulation, 86(7), 1281-1294.

*4:Walker, A. J. (1977). An efficient method for generating discrete random variables with general distributions. ACM Transactions on Mathematical Software (TOMS), 3(3), 253-256.