Who we are?
こんにちは。私達は justInCase という保険会社のスタートアップです。エンジニアはフラットな組織環境の中で密度の濃い毎日を過ごしており、「新しいことをやりたい！」、「自分の理想のチームを作り上げたい！」という意気込みがある方はふらっとオフィスに遊びに来て下さい。表には出せないけど面白いことをいろいろやってるのでお話しましょう！
さて、Rです。チームの中では少数民族で日常的な使用者は私しかおりません。前処理および探索的分析 with tidyverse な使い方が中心であり、プロダクションに組み込む実装ではほぼ python という使い分けになっています。そんな日常の中で R に潜むふとした疑問であった NA を今回取り上げます。なお会社設立時に勢いで書いた R のコードがオーパーツのような存在になりつつあるので、 R を使って前処理と前処理や前処理などを鍛えたい方は、ぜひぜひいらして下さい。
要約
- R の
NAは言語レベルで論理型、整数型、数値型、複素数型、文字列型それぞれにNA,NA_integer_,NA_real_,NA_complex_,NA_character_が存在する。 - これらの実装を見るにはC言語の実装を見れば良い。
.Internal(inspect(,,,))と github 上の Cの実装を活用しよう。 - R の
NA_real_は C言語のNaNだけど下位bitに 1954 が仕込まれてるよ。
Before starting
さて、R は動的言語であり全てがオブジェクトという性質を持ちますが、オブジェクトの種類についてはここでは解説しません。R Internals をご参照下さい。R のC言語での実装は cons cell を CAR と CDR で操作のS式だらけなので、実は Lisper の皆様にとても相性が良いと言えませす（言えません）。
第一部 NAとはNAんぞや
まずは NA を見てみましょう。
> .Internal(inspect(NA))
@7fcef4b1cd08 10 LGLSXP g0c1 [NAM(3)] (len=1, tl=0) -2147483648
最初の @7fcef4b1cd08 はオブジェクトのアドレス。 10 は次の LGLSXP に対応する番号です。上の R Internals にもこの番号がありましたね。さて、これらの意味をどうやってそれを知ることができるのかって？ それは R のソースコードが github 上 に公開されているので、検索することでわかります。primitive な function は接頭語に do_ を付けた名前が C言語にあるので、その変数名を検索すれば探すことができます。今回は少し特殊で .Internal の第１引数に渡す call する関数に do_ を付ける、したがって do_inspect を検索することで中身を知ることができます。中身は inspect_treeのこの辺りです。
ちなみに .Internal(inspect(x,deep=-1,pvec=5)) と第２引数はオブジェクトの深さ, 第３引数は配列の表示の長さを指定できます。 deep=−１ は無制限再帰の指定なので、どでかいオブジェクト見る場合は deep=0 から少しづつはじめると良いでしょう。
Rprintf("@%lx %02d %s g%dc%d [", (long) v, TYPEOF(v), typename(v),
v->sxpinfo.gcgen, v->sxpinfo.gccls);
g0g1 は garbage collector (GC) の世代ですね。 GCには参照カウント、マーク＆スイープなどがありますが闇深くなるので言及しません。次の [NAM(3)] は NAMED(v) に対応しています。
if (NAMED(v)) { if (a) Rprintf(","); Rprintf("NAM(%d)",NAMED(v)); a = 1; }
NAMED とはこれまた 検索 すると #define NAMED(x) ((x)->sxpinfo.named) とあるので sxpinfo.named を調べるために struct sxpinfo_struct を見てみると。。。と、この調子だと終わらないので各自ご自分でお調べ下さい。なお NAMED はオブジェクトをコピーした際に変更の必要の有無をこのflagで判断することによりメモリーを節約する仕組みです（厳密には reference counting が導入されたという歴史がありますが、 # define REFCNT(x) ((x)->sxpinfo.named), #define NAMED(x) ((x)->sxpinfo.named) と本質的にこの sxpinfo.named の部分は使われていることは同じです。） https://developer.r-project.org/Refcnt.html https://cran.r-project.org/doc/manuals/r-release/R-ints.html#Rest-of-header
本題の NA に戻ります。
さて、肝心の -2147483648 に対応する部分は (int) LOGICAL_ELT(v, i) なので 論理型配列のi番目を (int) に型変換したものです。 inspect の対応するコード部分 で pvec で配列の表示の長さを制御していることが確認できます。
switch (TYPEOF(v)) { /* for native vectors print the first elements in-line */
case LGLSXP:
if (XLENGTH(v) > 0) {
unsigned int i = 0;
while (i < XLENGTH(v) && i < pvec) {
Rprintf("%s%d", (i > 0) ? "," : " ",
(int) LOGICAL_ELT(v, i));
i++;
}
if (i < XLENGTH(v)) Rprintf(",...");
}
break;
プログラミング歴長い人は -2147483648 見てすぐ勘付きますね。 -2^31 の整数値は R では 論理型 NA に予約されているので使用することができません。
> -2^31
[1] -2147483648
> -2147483648L
[1] -2147483648
警告メッセージ:
非整数値 2147483648L に L がついています; 実数値を使用します
その実態は arithmetic.c(L161-175) で定義されています
void attribute_hidden InitArithmetic()
{
R_NaInt = INT_MIN;
R_NaReal = R_ValueOfNA();
// we assume C99, so
#ifndef OLD
R_NaN = NAN;
R_PosInf = INFINITY;
R_NegInf = -INFINITY;
#else
R_NaN = 0.0/R_Zero_Hack;
R_PosInf = 1.0/R_Zero_Hack;
R_NegInf = -1.0/R_Zero_Hack;
#endif
}
はい。 INT_MIN が確認できました。まずは論理型の NA の実体をつかむプロセスが身についたと思います。
さて、このペースでやっていと、チームの人々から仕事しろと怒らてしまうので NA_real_ の解説は続けますが、それ以外は皆さんの宿題にしましょう（よくある丸投げ）。
> .Internal(inspect(NA_real_))
@7fcef4495430 14 REALSXP g0c1 [NAM(3)] (len=1, tl=0) nan
> .Internal(inspect(NA_complex_))
@7fcef234ad08 15 CPLXSXP g0c2 [NAM(3)] (len=1, tl=0)
> .Internal(inspect(NA_character_))
@7fcef32d8820 16 STRSXP g0c1 [NAM(3)] (len=1, tl=0)
@7fcef2807200 09 CHARSXP g0c1 [MARK,gp=0x20] [cached] "NA"
# 参考
> .Internal(inspect("NA"))
@7fcef32d8778 16 STRSXP g0c1 [NAM(3)] (len=1, tl=0)
@7fcef4a58080 09 CHARSXP g0c1 [gp=0x60] [ASCII] [cached] "NA"
おっと衝撃です。 NA_real_ はCの nan なのか？逆に Rの NaN を調べてみましょう。
> .Internal(inspect(NaN))
@7fcef32d8660 14 REALSXP g0c1 [NAM(3)] (len=1, tl=0) nan
あれ、これも同じ nan ですね。なんなんでしょうか？（これが言いたかっただけの記事です）。 identica() でチェックするも当然異なりますね。てことはソースコードへの旅にでかけましょう。
> identical(NaN, NA_real_)
[1] FALSE
第２部 nan は nan なのか
さきほどの arithmetic.c をもう一度見てみましょう。
R_NaReal = R_ValueOfNA();
// we assume C99, so
#ifndef OLD
R_NaN = NAN;
この R_ValueOfNA() 探すとまたもや arithmetic.c に出てきました！
https://github.com/wch/r-source/blob/5a156a0865362bb8381dcd69ac335f5174a4f60c/src/main/arithmetic.c#L108-L116
static double R_ValueOfNA(void)
{
/* The gcc shipping with Fedora 9 gets this wrong without
* the volatile declaration. Thanks to Marc Schwartz. */
volatile ieee_double x;
x.word[hw] = 0x7ff00000;
x.word[lw] = 1954;
return x.value;
}
そもそも NaN は IEEE754 で定義されていますが、浮動小数点表現では wikipedia によると
a bit-wise IEEE floating-point standard single precision (32-bit) NaN would be:
s111 1111 1xxx xxxx xxxx xxxx xxxx xxxx where s is the sign (most often ignored
in applications) and the x sequence represents a non-zero number (the value zero
encodes infinities)
とあります。もうすこし分かりやすく日本語のwikipediから IEEE754 の32bit 浮動小数点の図を借りて確認すると、sign は無視、exponent (8bits) 部分が全て1、 fraction (23bits) 部分は non-zero が NaN と定義されています。
(出典 wikipedia: https://upload.wikimedia.org/wikipedia/commons/thumb/d/d2/Float_example.svg/1180px-Float_example.svg.png)
さて ieee_double は 同ファイルの L83-87 に
typedef union
{
double value;
unsigned int word[2];
} ieee_double;
と union で定義され、 R_ValueOfNA() では unsigned int word [2] を使いなので上位と下位に 0x7ff00000 と 1954 をぶっこんだ NaN ですね。（この hw, lw はもちろん endianness で判定してます）。これが Rの NA_real_ すなわち C言語のNaNの下位の任意値に 1954 が入っているものです。bit 表現には library(R.utils::intToBin) を使うことにより書くにできます。当然 IEEE754 の NaN の定義を満たしている表現です。
> library(R.utils)
> intToBin(0x7ff00000 + 1954)
[1] "1111111111100000000011110100010"
この 1954 については 2004年の R mailing list に言及があります。 https://stat.ethz.ch/pipermail/r-help/2004-September/057043.html こういう easter egg を入れられるのが言語設計者の魅力でもありますね。Rマウンティング取りたい方、いつかこの知識を披露して下さい。
要約（再掲）
- R の
NAは言語レベルで論理型、整数型、数値型、複素数型、文字列型それぞれにNA,NA_integer_,NA_real_,NA_complex_,NA_character_が存在する。 - これらの実装を見るにはC言語の実装を見れば良い。
.Internal(inspect(,,,))と github 上の Cの実装を活用しよう。 - R の
NA_real_は C言語のNaNだけど下位bitに 1954 が仕込まれてるよ。
最後に
こういうマニアックな話が大好きな方、ぜひ justInCase へ遊びに来て下さい。https://www.wantedly.com/companies/justincase