NAとはNAんぞや (R)

## Who we are?
こんにちは。私達は justInCase という保険会社のスタートアップです。エンジニアはフラットな組織環境の中で密度の濃い毎日を過ごしており、「新しいことをやりたい！」、「自分の理想のチームを作り上げたい！」という意気込みがある方はふらっとオフィスに遊びに来て下さい。表には出せないけど面白いことをいろいろやってるのでお話しましょう！

さて、Rです。チームの中では少数民族で日常的な使用者は私しかおりません。前処理および探索的分析 with tidyverse な使い方が中心であり、プロダクションに組み込む実装ではほぼ python という使い分けになっています。そんな日常の中で R に潜むふとした疑問であった `NA` を今回取り上げます。なお会社設立時に勢いで書いた R のコードがオーパーツのような存在になりつつあるので、 R を使って前処理と前処理や前処理などを鍛えたい方は、ぜひぜひいらして下さい。

## 要約

1. R の `NA` は言語レベルで論理型、整数型、数値型、複素数型、文字列型それぞれに `NA`, `NA_integer_`, `NA_real_`, `NA_complex_`, `NA_character_` が存在する。
1. これらの実装を見るにはC言語の実装を見れば良い。 `.Internal(inspect(,,,))` と github 上の Cの実装を活用しよう。
1. R の `NA_real_` は C言語の `NaN` だけど下位bitに 1954 が仕込まれてるよ。

## Before starting
さて、R は動的言語であり全てがオブジェクトという性質を持ちますが、オブジェクトの種類についてはここでは解説しません。[R Internals](https://cran.r-project.org/doc/manuals/r-release/R-ints.html#SEXPTYPEs) をご参照下さい。R のC言語での実装は cons cell を `CAR` と `CDR` で操作のS式だらけなので、実は Lisper の皆様にとても相性が良いと言えませす（言えません）。

## 第一部 NAとはNAんぞや
まずは `NA` を見てみましょう。

```R
> .Internal(inspect(NA))
@7fcef4b1cd08 10 LGLSXP g0c1 [NAM(3)] (len=1, tl=0) -2147483648
```

最初の `@7fcef4b1cd08` はオブジェクトのアドレス。 `10` は次の `LGLSXP` に対応する番号です。上の R Internals にもこの番号がありましたね。さて、これらの意味をどうやってそれを知ることができるのかって？ それは [R のソースコードが github 上](https://github.com/wch/r-source/) に公開されているので、検索することでわかります。primitive な function は接頭語に `do_` を付けた名前が C言語にあるので、その変数名を検索すれば探すことができます。今回は少し特殊で `.Internal` の第１引数に渡す call する関数に `do_` を付ける、したがって `do_inspect` を検索することで中身を知ることができます。中身は [inspect_tree](https://github.com/wch/r-source/blob/5f0affa2c7016e054f3eb4b64e247d428a6477dd/src/main/inspect.c#L96)のこの辺りです。
ちなみに `.Internal(inspect(x,deep=-1,pvec=5))` と第２引数はオブジェクトの深さ, 第３引数は配列の表示の長さを指定できます。 `deep=−１` は無制限再帰の指定なので、どでかいオブジェクト見る場合は `deep=0` から少しづつはじめると良いでしょう。

```C
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)` に対応しています。

```R
if (NAMED(v)) { if (a) Rprintf(","); Rprintf("NAM(%d)",NAMED(v)); a = 1; }
```

`NAMED` とはこれまた [検索](https://github.com/wch/r-source/blob/3b71bbc1f2f1f6aca1158ec47f1735daedf9c67b/src/include/Rinternals.h#L315) すると `#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` の[対応するコード部分](https://github.com/wch/r-source/blob/5f0affa2c7016e054f3eb4b64e247d428a6477dd/src/main/inspect.c#L150-L160) で `pvec` で配列の表示の長さを制御していることが確認できます。

```R
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 に予約されているので使用することができません。

```R
> -2^31
[1] -2147483648
> -2147483648L
[1] -2147483648
 警告メッセージ: 
 非整数値 2147483648L に L がついています; 実数値を使用します 
```

その実態は [arithmetic.c(L161-175)](https://github.com/wch/r-source/blob/5a156a0865362bb8381dcd69ac335f5174a4f60c/src/main/arithmetic.c#L163) で定義されています

```C
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_` の解説は続けますが、それ以外は皆さんの宿題にしましょう（よくある丸投げ）。

```R
> .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` を調べてみましょう。

```R
> .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 をもう一度見てみましょう。

```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

```C
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 によると

```text
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` と定義されています。
![test](https://upload.wikimedia.org/wikipedia/commons/thumb/d/d2/Float_example.svg/1180px-Float_example.svg.png)
(出典 wikipedia: https://upload.wikimedia.org/wikipedia/commons/thumb/d/d2/Float_example.svg/1180px-Float_example.svg.png)

さて `ieee_double` は 同ファイルの L83-87 に

```C
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 の定義を満たしている表現です。

```R
> 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マウンティング取りたい方、いつかこの知識を披露して下さい。




## 要約（再掲）

1. R の `NA` は言語レベルで論理型、整数型、数値型、複素数型、文字列型それぞれに `NA`, `NA_integer_`, `NA_real_`, `NA_complex_`, `NA_character_` が存在する。
1. これらの実装を見るにはC言語の実装を見れば良い。 `.Internal(inspect(,,,))` と github 上の Cの実装を活用しよう。
1. R の `NA_real_` は C言語の `NaN` だけど下位bitに 1954 が仕込まれてるよ。


## 最後に

こういうマニアックな話が大好きな方、ぜひ justInCase へ遊びに来て下さい。https://www.wantedly.com/companies/justincase