matlab: なんちゃってcomplex

matlab (今回は 6.5.1 でしか試してません)は複素数について直感に反する挙動を示すことがある．例を見てみよう．

>> A = repmat(complex(0),[10,1]);
>> whos A
  Name      Size                   Bytes  Class
  A        10x1                      160  double array (complex)
Grand total is 10 elements using 160 bytes
>> A(end) = complex(1,0);
>> whos A
  Name      Size                   Bytes  Class
  A        10x1                       80  double array
Grand total is 10 elements using 80 bytes

complex じゃなくなってる．およ？一応確認すると，A に代入した値は実数でなく複素数と認識されているはずだ．

>> clear a; a=complex(1,0); whos a
  Name      Size                   Bytes  Class
  a         1x1                       16  double array (complex)
Grand total is 1 element using 16 bytes

となるからだ．ただし，入力引数に対する inputname のように，出力引数についての情報取得機構が関数 complex 内で働いているとすると，よくわからない話になるかもしれない．

一つ考えられる仮説は，「虚部が全て 0 の配列は実数型にされてしまう」というものだ．実際，その方がメモリ効率やら計算効率やらよくなる気はする．

念のため，虚部が 0 でない複素数を代入すれば，配列は実数型に変換されることは無い．

>> A = repmat(complex(0),[10,1]);
>> A(end)=complex(0,1);
>> whos A
  Name      Size                   Bytes  Class
  A        10x1                      160  double array (complex)
Grand total is 10 elements using 160 bytes

そして，要素を虚部が 0 の複素数に順に置き換えていくと，

>> A=complex([1,2],[3,4]);
>> A(1)=complex(real(A(1))); whos A
  Name      Size                   Bytes  Class
  A         1x2                       32  double array (complex)
Grand total is 2 elements using 32 bytes
>> A(2)=complex(real(A(2))); whos A
  Name      Size                   Bytes  Class
  A         1x2                       16  double array
Grand total is 2 elements using 16 bytes

となった．ここまでは前述の仮説で当たっていそうである．

しかし，もちろん話はこれで終わらない．次の例を見てみよう．

>> A=zeros(1,2);
>> A(:)=complex(0);
>> whos A
  Name      Size                   Bytes  Class
  A         1x2                       32  double array (complex)
Grand total is 2 elements using 32 bytes

なんと配列は複素数のままであり，これは前述の仮説と矛盾する．ところがさらに，代入の方法を変えると，

>> A=zeros(1,2);
>> A(1:end)=complex(0);
>> whos A
  Name      Size                   Bytes  Class
  A         1x2                       16  double array
Grand total is 2 elements using 16 bytes

となって前述の仮説が復活する．どうなってるのという話である．

ここで新たに，｢全代入(whole substition とか？)｣という概念を考えてみる．
全代入と考えるのは次の４つ．

1. A = B;
2. A(:) = B;
3. A = B(:);
4. A(:) = B(:);
   

変数間の代入で虚部 0 でも complex を維持するには，どうやら厳密にこの４形式しか許されないようだ．
たとえば２次元配列で A(:,:) としてもダメだったし，列ベクトルで A(:,1) としてもダメだった．
この概念を利用して立てた仮説がこれだ：
「全代入では愚直な代入が行われ，それ以外では値の取り出し・代入いずれかで可能なら虚部が 0 の複素数は実数に落ちる」
いまのところこの仮説で説明できるような気がする．
この他に，左辺が全参照（全代入と同じ気分）で右辺が「適当な」式なら，やはり実数に落ちないようだ．以下の例を参照のこと．

>> A=complex(rand(3));
>> A(:)=complex(0); whos A
  Name      Size                   Bytes  Class
  A         3x3                      144  double array (complex)
Grand total is 9 elements using 144 bytes
>> A(:)=complex(zeros(3)); whos A
  Name      Size                   Bytes  Class
  A         3x3                      144  double array (complex)
Grand total is 9 elements using 144 bytes

いや～．重箱の隅つつくといろいろ出てくるなぁ．
ちなみに，大方の予想通り

>> A=complex(rand(3));
>> A(:,:)=complex(0); whos A

  Name      Size                   Bytes  Class
  A         3x3                       72  double array
Grand total is 9 elements using 72 bytes
>> A=complex(rand(3));
>> A(:)=complex(zeros(3)) + 1; whos A
  Name      Size                   Bytes  Class
  A         3x3                       72  double array
Grand total is 9 elements using 72 bytes

となる．もういいって？
今回の現象は，subsref, subsasgn, subsindex と型キャストの実装のコンボで生じるんだろうけど，意図された仕様なんですかねぇ．