nptclのブログ

Common Lisp処理系nptの開発メモです。https://github.com/nptcl/npt

Common LispとC言語で虚数部のゼロ符号に違いが出る

ずっと悩んでいたことがわかったのでメモします。
例えば下記の関数

log(1-z)

のzに10を入れると、Common LispC言語では結果が違います。

Common Lispの例

(defun log1-z (z)
  (log (- 1.0 z)))
(log1-z 10)
  -> #C(2.1972246 3.1415927)

C言語の例

#include <stdio.h>
#include <complex.h>
#include <math.h>

complex double log1_z(complex double z)
{
    return clog(1.0 - z);
}

int main()
{
    complex double z = log1_z(10.0);
    printf("%f, %f\n", creal(z), cimag(z));
    return 0;
}
$ cc main.c -lm
$ ./a.out
2.197225, -3.141593

結果

Common Lisp
2.1972246  3.1415927
C言語
2.197225  -3.141593

虚数部の符号が違いますね。
それだけなんですが、なぜ?

確認していったら、どうも関数logに渡される引数の 虚数部ゼロの符号が違っているからのようです。 下記の式の部分です。

1.0-z

Common Lispでは

1.0 - 10.0 = -9.0

と当たり前の数を出しますが、 C言語では引数がそもそも複素数で指定されているため、

1.0 - 10.0 = -9.0 -0.0*I

のようになります。 このゼロにマイナスがあるかないかで、logの値が変わっていました。 数学に全然詳しくないんですが、これはどっちも正しいでいいんですかね。 ダメだとしても、差異はどうしようもなさそうです。

言語レベルでの話をしますと、Common Lispでは簡単にはC言語の様な挙動に できないと思います。
関数呼び出し時に

(log1-z 10)

(log1-z #c(10 0))

としたところで、

#c(10 0)

10

に変換されてから関数に渡されてしまいます。

【追記】間違いです。 #c(10 0)ではなく浮動小数#c(10.0 0.0)とすればcomplex型のままになります。 以降の内容にはそれほど影響がなかったため、内容はそのままにします。

ではlog1-zの関数定義を下記のように変更したらどうか。

(defun log1-z (z)
  (declare (type complex z))
  (log (- 1.0 z)))

ダメですね。 declare typeは型チェックを行う役目しかありません。 sbcl, clisp, cclで確認しましたが期待通りの動きはしません。

ただし、負ゼロの値を渡すことはできる場合はあります。 下記をsbclとcclにて確認しました。

* (log -10.0)
#C(2.3025851 3.1415927)
* (log #c(-10.0 -0.0))
#C(2.3025851 -3.1415927)

ただしclispではダメ。

まあ、本当に必要なら符号見て場合わけすりゃいい話なんですが、 標準関数内に上記のコードが埋め込まれている場合があります。 上記の差異はarcsinを確認しているときに見つかりました。

Common Lispでのasin(30)
  #C(1.5707964 -4.0940666)

C言語でのcasin(30.0)
  1.570796, 4.094067

arcsinの定義を追っていくと、

arcsin(z) = -i*log(i*z + sqrt(1-z*z)) 

となり、上記sqrt部が

sqrt(1-z*z) = exp(log(1-z*z) / 2.0)

となるため、log(1-z*z)が含まれるので、虚数部ゼロの符号に差異が生じるというわけです。

何が正しいのか全くわからん。