コグノスケ

三角関数のfloat版

目次: C言語とlibc

標準Cライブラリにはdoubleを返す三角関数（sin(), cos(), tan()）とfloatを返す三角関数（sinf(), cosf(), tanf()）が定義されています。

標準Cライブラリの一つの実装であるmuslのコードを見ると、sinf, cosf, tanfの計算にdouble演算を内部で使っています。これは基になったFreeBSDのlibmと同じ実装です。PCはfloatでもdoubleでも関係なく速いんですが、doubleをハードで扱えない貧相なプロセッサには優しくない作りです。

もう一つの実装であるnewlibのコードを見ると、double版のsin, cosこそFreeBSDの実装と同じですが、float版のsinf, cosfはfloatを使ったコードが独自に追加されていて、貧相なプロセッサにも優しい作りになっています。組み込みにやたら使われる実績は伊達じゃないですね。

じゃあmuslにnewlibのfloat版のsinf, cosfを移植すれば、doubleが苦手なプロセッサでも速くなるのでは？と思いました。

テストを先に書こう

コードを触る前に、それぞれの実装の素性を調べておこうと思います。テスト方法は、

• 期待値: glibcのsin, cos, tan（double版）をfloatに変換した結果
• 判定方法: sinf, cosfは1の誤差を許す、tanfは3の誤差を許す
• 比較範囲: floatを32bitとして、全値域（＝約43億パターン）

どうしてtanfだけ判定が甘いかというと、正しい値がわからなかったからです。なぜかglibcの実装も誤差1に収まっていません。どういうことなの……。

テストのコードはGitHubに置き（リンク）ました。特に難しい点はありませんが、muslとnewlibから三角関数を拝借するところは、やや面倒かもしれません。

現在のプロセッサは超速いし、問題の性質上マルチスレッド化も簡単ですから、32並列くらいで頑張れば32bit全域を調査しても3分もかかりません。楽勝ですね〜。

最初にテストして良かった

テスト結果は、当然、全て一致かと思いきや、そんなことはなかった。最初にテストしておいて良かったですね。

良い方から言うとmuslは内部でdoubleで演算しているからか、結果もパーフェクトでした。

一方のnewlibはcosfだけ変な値を返します。32bit全域を試してわずか6パターンです。

cos,cosf_newlib: NG : x:3fc90fe0 f:1.570797 d:1, exp:b52bbbd3 -0.000001, res:b52bbbd0 -0.000001
cos,cosf_newlib: NG : x:3fc90fe1 f:1.570797 d:1, exp:b54bbbd3 -0.000001, res:b54bbbd0 -0.000001
cos,cosf_newlib: NG : x:3fc90fe2 f:1.570797 d:1, exp:b56bbbd3 -0.000001, res:b56bbbd0 -0.000001

cos,cosf_newlib: NG : x:bfc90fe0 f:-1.570797 d:1, exp:b52bbbd3 -0.000001, res:b52bbbd0 -0.000001
cos,cosf_newlib: NG : x:bfc90fe1 f:-1.570797 d:1, exp:b54bbbd3 -0.000001, res:b54bbbd0 -0.000001
cos,cosf_newlib: NG : x:bfc90fe2 f:-1.570797 d:1, exp:b56bbbd3 -0.000001, res:b56bbbd0 -0.000001

正負を考慮（浮動小数点は最上位ビットが符号を示すビット）すると、実質3パターンで変な値が返ることがわかります。

• 正: 0x3fc90fe0, 0x3fc90fe1, 0x3fc90fe2
• 負: 0xbfc90fe0, 0xbfc90fe1, 0xbfc90fe2

誤差は3でした。ほぼ合ってます、おしい。誤差が出ることも不思議ですが、sinfは合っていてcosfだけ値がズレるのも不思議です。

GCCを調べる - その8-3 - レジスタconstraint判定

目次: GCC

インラインアセンブラで "v" constraintsを指定すると、何も実装していない場合はimpossible constraint in 'asm' と怒られました。レジスタのconstraintsだけ足すとinconsistent operand constraints in an asmと怒られるはずです。エラーをチェックしている箇所は、

static bool
curr_insn_transform (bool check_only_p)
{

...

  if (process_alt_operands (reused_alternative_num))  //★★これが成立してalt_p = trueが期待値だが
    alt_p = true;

...

  if (! alt_p && ! sec_mem_p)
    {
      /* No alternative works with reloads??  */
      if (INSN_CODE (curr_insn) >= 0)
	fatal_insn ("unable to generate reloads for:", curr_insn);
      error_for_asm (curr_insn,
		     "inconsistent operand constraints in an %<asm%>");  //★★ここに到達しエラーが出る
      lra_asm_error_p = true;
      /* Avoid further trouble with this insn.  Don't generate use
	 pattern here as we could use the insn SP offset.  */
      lra_set_insn_deleted (curr_insn);
      return true;
    }

...

このcurr_insn_transform() 関数はやたら長くて（700行）訳のわからない構造です。うまく行く場合（rなどを渡したとき）を観察すると、alt_pがtrueになるのが期待値と思われます。幸いなことにalt_pの設定は一箇所だけ、条件もprocess_alt_operands() 関数だけです。

そう思ってprocess_alt_operands() 関数を見ると、これがまたもの凄い実装で、目を覆いたくなります（1200行！！）。GCC見ていると、クソコードには事欠かないです。これはひどい。

コードの一部を抜粋しても全く意味不明で、そもそもこの関数自体がかなりゴチャゴチャで意味不明です。全て追うのは不可能です。なので"r" がどの辺りを通るかをもって、当たりを付けました。下記のところが分岐点になっているようです。

static bool
process_alt_operands (int only_alternative)
{

...

	  do
	    {
              //★★pは "=&v" が入っていて、cに先頭から一文字ずつ取って解析している
	      switch ((c = *p, len = CONSTRAINT_LEN (c, p)), c)
		{
		case '\0':
		  len = 0;
		  break;

...

		default:
		  cn = lookup_constraint (p);  //★★ 'v' に対しては、CONSTRAINT_vが返る
		  switch (get_constraint_type (cn))
		    {
		    case CT_REGISTER:
		      cl = reg_class_for_constraint (cn);  //★★CONSTRAINT_vに対してはVP_REGSが返る
		      if (cl != NO_REGS)
			goto reg;    //★★このジャンプで飛ぶ
		      break;

...

		reg:
		  if (mode == BLKmode)
		    break;
		  this_alternative = reg_class_subunion[this_alternative][cl];
		  this_alternative_set |= reg_class_contents[cl];  //★★どこかでみたreg_class_contentsが登場
		  if (costly_p)
		    {
		      this_costly_alternative
			= reg_class_subunion[this_costly_alternative][cl];
		      this_costly_alternative_set |= reg_class_contents[cl];
		    }
		  winreg = true;
		  if (REG_P (op))
		    {
		      if (hard_regno[nop] >= 0
			  && in_hard_reg_set_p (this_alternative_set,
						mode, hard_regno[nop]))  //★★これが成立しない
			win = true;  //★★少なくともwin = trueにならないと関数が失敗を返す（条件は他にもあるが）
		      else if (hard_regno[nop] < 0
			       && in_class_p (op, this_alternative, NULL))
			win = true;
		    }
		  break;
		}

...

	    }
	  while ((p += len), c);  //★★基本は次の文字に行くが、スキップすることもある模様

どこかでみたアイツです。このエラーはreg_class_contentsを見に行った結末に起きているようです。

試してみたら、色々おかしい

REG_CLASS_CONTENTSを正しく設定すると、下記のコードがコンパイルできるはずです。雰囲気を出すためRISC-Vのベクトル命令を書いていますが、ぶっちゃけコンパイラは命令を全く見ないので、実はabcdでも何でも通ります。コンパイルのみ（*.sを出力）であればアセンブラすら要りません（※）。

// a.c

void _start()
{
	int b[100];
	int v;

	__asm__ volatile ("vlw.v %0, %1\n"
		: "=&v"(v) : "A"(b[10]));
}

ビルドして、逆アセンブルしてみます。

$ riscv32-unknown-elf-gcc -Wall -g -march=rv32gcv -mabi=ilp32f -nostdlib -O2 a.c

$ riscv32-unknown-elf-objdump -drS a.out

a.out:     file format elf32-littleriscv

Disassembly of section .text:

00010054 <_start>:
void _start()
{
   10054:       7165                    addi    sp,sp,-400
        int b[100];
        int v;

        __asm__ volatile ("vlw.v %0, %1\n"
   10056:       103c                    addi    a5,sp,40
   10058:       1207e007                vlw.v   v0,(a5)
                : "=&v"(v) : "A"(b[10]));
}
   1005c:       6159                    addi    sp,sp,400
   1005e:       8082                    ret

それらしきベクトルレジスタ（v0）が出力されているようです。めでたし、めでたし。と言いたいところですが、実は全然ダメです。

• 変数がintなのでsizeof(v) が4になる、ベクトルを扱いたい
• 最適化オプションをO0にするとコンパイラがinternal errorを出す

まだまだ改善の余地があります。これも今後、調べていこうと思います。

（※）もしアセンブルまで実行したければ、RISC-VのGitHubにあるbinutilsを使ってください（GitHubへのリンク）。ビルド方法はUpstreamのコードとほぼ同じ（2019年4月19日の日記参照）です。唯一の違いはconfigure時に --with-system-readlineを付けないと、readlineがないと言われてエラーになる点です。