はじめに

「integer64で定義されたベクトルをmatrixに変換したら値がおかしくなった」という話を聞いたので原因を調査してみました。あえて結論を最初に書かないので、忙しい人はまとめを参照してください。

再現可能な例

library(bit64)

int64 <- as.integer64(c(1, 2, 3))
int64

integer64
[1] 1 2 3

int64 |>
  as.matrix()

              [,1]
[1,] 4.940656e-324
[2,] 9.881313e-324
[3,] 1.482197e-323

わけがわからないくらい小さな数字に変化してしまっています。暇な人はスクロールを止めて原因を考えてみましょう。

integer64ってなんだ

なんとなくinteger64が怪しそうだし、numericに変換すれば直るんじゃね？ｗくらいのノリでコメントしたら予想通り解決しました。

int64 |>
  as.numeric() |>
  as.matrix()

     [,1]
[1,]    1
[2,]    2
[3,]    3

原因調査

検索すると似たような問題で躓いている人が他にも結構いるみたいです¹。動けばOK！ｗを卒業したいという気持ちもあるので原因を考えていきます。

numericに変換すると期待通りに動作する理由を考える

まずはas.numericが何をやっているかを見ていきます。

as.numeric

function (x, ...)  .Primitive("as.double")

as.numericはas.doubleを呼んでいるんですね。bit64を読み込んだときにas.doubleがマスクされていたのでメソッドを確認してみます。

getS3method("as.double", "integer64")

function (x, keep.names = FALSE, ...) 
{
    ret <- .Call(C_as_double_integer64, x, double(length(x)))
    if (keep.names) 
        names(ret) <- names(x)
    ret
}
<bytecode: 0x00000128027d9408>
<environment: namespace:bit64>

.Callって何？という気持ちになります。ドキュメントを見ると内部関数を読み込むものらしいです。bit64のソースコードがGitHubに置いてあったので²検索を掛けてみるとそれっぽい関数がヒットしました。

SEXP as_double_integer64(SEXP x_, SEXP ret_){
  long long i, n = LENGTH(x_);
  long long * x = (long long *) REAL(x_); 
  double * ret = REAL(ret_);
  double rmax = pow(FLT_RADIX, DBL_MANT_DIG) - 1;
  double rmin = -rmax;
  Rboolean naflag = FALSE;
  for (i=0; i<n; i++){
    if (x[i]==NA_INTEGER64)
      ret[i] = NA_REAL;
      else{
        if (x[i]<rmin || x[i]>rmax)
            naflag = TRUE;
        ret[i] = (double) x[i];
      }
  }
  if (naflag)warning(INTEGER64_TODOUBLE_WARNING);
  return ret_;
}

本質的にはret[i] = (double) x[i];の行でlong long型をdouble型にキャストしているだけっぽいですね。これでうまくいく理由はわかりました。次は上手くいかない方の理由を考えていきます。

matrixはなにをしているか

matrixのメソッドから確認してみましょう。

methods("matrix")

no methods found

なかったです。次に中身を見てみることにします。

matrix

function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL) 
{
    if (is.object(data) || !is.atomic(data)) 
        data <- as.vector(data)
    .Internal(matrix(data, nrow, ncol, byrow, dimnames, missing(nrow), 
        missing(ncol)))
}
<bytecode: 0x000001287114f568>
<environment: namespace:base>

またC言語の予感。ぐえ～。

同様に実装を検索してみたところ、みつかりましたが結構長いです³。型まわりが怪しいはずなのでまずはこの部分を読んでみます。

switch(TYPEOF(vals)) {
    case LGLSXP:
    case INTSXP:
    case REALSXP:
    case CPLXSXP:
    case STRSXP:
    case RAWSXP:
    case EXPRSXP:
    case VECSXP:
        break;
    default:
        error(_("'data' must be of a vector type, was '%s'"),
            R_typeToChar(TYPEOF(vals)));
}

TYPEOF()⁴でmatrix()の引数の型をチェックしていそうです。また、行列を作成するためにTYPEOF(vals)とおなじ型でメモリを確保しています。

PROTECT(ans = allocMatrix(TYPEOF(vals), nr, nc));

そういえばRにもtypeof()が実装されていた気がします。確認してみましょう。

typeof

function (x) 
.Internal(typeof(x))
<bytecode: 0x0000012873086568>
<environment: namespace:base>

matrix()で呼んでいる内部関数と同じっぽいですね。integer64型のtypeを確認してみます。

int64 |>
  typeof()

[1] "double"

おやおや。double型が出てきました。どうもmatrixはinteger64をdoubele型として解釈しているようです。こいつが悪さしていそうに思えます。次はinteger64をdoubele型として解釈するとどうなるのかを考えていくことにします。

integer64型の値をdoubel型として解釈してみる

C言語の授業を思い出す時間がやってきました。integer64、いわゆるlong long intは先頭に1bitの符号ビット、残りの63bitが値を表す部分です。整数の1を内部表現で表すと次のようになります。

00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000001

この内部表現をdouble型として解釈するとどうなるか考えます。double型は1bitの符号ビット、11bitの指数部、52bitの仮数部で構成されます。今回のケースだと次の通りです。

符号ビット: 0 (+) 指数部: 00000000000 (？) 仮数部: 0000000000000000000000000000000000000000000000000001

指数部が全て0になることなんてあったっけ？となってしらべてみたら、ありました（無知）。非正規化数と呼ぶらしいです⁵。 Wikipediaに記載された定義通り計算してみます。

\[ \begin{align*} 値 &= 符号 \cdot 2^{指数部} \cdot 仮数部 \\ &= (+1) \cdot 2^{-1022} \cdot (1.0 \cdot 2^{-52}) \\ &= 2^{-1022} \cdot 2^{-52} \\ &= 2^{-1074} \end{align*} \]

2^(-1074)

[1] 4.940656e-324

勝ちです。

まとめ

integer64型はdouble型で実装されている
bit64の読み込み時に多くの関数はinteger64型に対応するメソッドが追加される
matrixはメソッドがなく、引数の型をもとにしてRの内部関数で処理を行う
結果として、matrixの内部でinteger64型をdouble型と誤認してしまうため大きく異なる値に変換されてしまう

matrixに限った問題ではないと思います。integer64型を利用するときは気を付けましょう！

bit64パッケージのinteger64型には注意しよう