warnの出力抑制はどうやって実現しているのか
warn
の出力を抑制する方法
(handler-bind ((warning #'muffle-warning)) (warn "Hello"))
やった!
warn
関数とその周辺の実装を考える
Common Lispのwarning
の話題です。
warn
関数では警告を出力できますが、
その出力を抑止したい場合はどうしたらいいでしょうか。
先行して回答を示したように、handler-bind
とmuffle-warning
を組み合わせて
設定することでwarn
の出力を抑止することができます。
つまり、単純に
(warn "Hello")
と実行すると、
WARNING: Hello
みたいな出力が出るのですが、
(handler-bind ((warning #'muffle-warning)) (warn "Hello"))
では何も出力されません。
これは一体どのような仕組みで実行されているのでしょうか?
handler-bind
とmuffle-warning
という組み合わせから、
condition
とrestart
が手を組んでいることがわかると思います。
たかがwarning
抑止だと舐めてかかれるようなモノじゃなく、
なんか知らんが結構複雑ですので、ちゃんと説明していきます。
muffle-warning
とは何か
muffle-warning
はrestart
関数と呼ばれています。
最後に使い方を説明しますが、warning
を中断させる機能があります。
restart
関数は他にもあり、abort
関数やらcontinue
関数がそれにあたります。
restart
関数の内容は、多少の差異はあるものの大体決まっており、
find-restart
とinvoke-restart
の2つにより構成されています。
muffle-warning
の実装例である、muffle-warning!
を下記に示します。
(defun muffle-warning! (&optional condition) (let ((restart (find-restart 'muffle-warning! condition))) (if restart (invoke-restart restart) (error (make-condition 'control-error)))))
たぶんどの処理系でもこんな感じに実装されていると思います。
find-restart
で探して、見つかったならそれをinvoke
するというもの。
つまりは、warning
というcondition
を受け取ったら、
muffle-warning
より、再起動関数が呼び出されることになります。
逆にいうなら、warning condition
をsignal
で呼び出す前には、
必ずrestart-bind
かrestart-case
にてmuffle-warning restart
を用意してあげなければ
control-error
になってしまいます。
この要件が、まずはwarn
関数を作成する際に必要になります。
注意点としては、restart
関数は自分の関数名と同じ名前のrestart
を探します。
今の場合、muffle-warning
関数は、muffle-warning
という名前のrestart
を
探してinvoke
しているのです。
関数名とrestart
の名前はわざと同じにしているのですが、
同じsymbol
でも無関係であることを覚えておいた方がいいと思います。
自作のwarning condition
を作る
テスト用にwarning condition
を作りましょう。
要件は1つ。format
の2要素である、format
文字列と引数を受け取れることです。
slot
を作るのが面倒なのでsimple-condition
を継承します。
名前はwarning!
とします。
(define-condition warning! (simple-condition) ()) (defmethod print-object ((instance warning!) stream) (let ((format (simple-condition-format-control instance)) (args (simple-condition-format-arguments instance))) (apply #'format stream format args)))
print-object method
は、例えばprinc
なんかで引数として与えられたときに、
どういう文字列を出力するかを指定するものです。
深く考えずに、warning!
が受け取った引数から、format
をそのまま出力してやります。
例えばこんな感じ。
(let ((inst (make-condition 'warning! :format-control "HELLO" :format-arguments nil))) (format t "~S~%" inst)) → HELLO
自作のwarn
関数を作る
では、上記の要件から、自作のwarn
関数である、warn!
関数を作成することを考えます。
要件は2つ。
muffle-warning! restart
を用意する。warning! condition
を実行する。
warn!
関数の実装は下記のようになります。
(defun warn! (format &rest args) (restart-case (signal (make-condition 'warning! :format-control format :format-arguments args)) (muffle-warning! ())))
それではwarn!
関数を次のように実行した場合はどうなるでしょうか。
(handler-bind ((warning! #'muffle-warning!)) (warn! "Hello"))
warn!
関数は、まずはmuffle-warning! restart
を用意したあと、
warn! instance
に出力内容を設定して、signal
でcondition
を呼び出します。
するとhandler-bind
により、muffle-warning!
関数が呼び出されます。
muffle-warning!
関数はinvokeによりwarn!
関数に戻り、restart-case
が実行されます。
restart-case
の本体には何も記載がないものの、
bind
ではなくcase
であるため、warn!
関数に制御がそのまま残り、
そして何もしないでwarn!
関数が終了します。
流れは次のような感じになります。
warn!
起動signal
によりhandler-bind
のwarning!
起動muffle-warning!
関数実行restart-case
のmuffle-warning!
起動restart-case
終了warn!
終了
これは本来のwarn
関数の挙動と同じになります。
それではhandler-bind
無しでwarn!
を呼び出した場合はどうなるでしょうか。
残念ながらこのままでは、warning! condition
のhandler
が存在しないため、
signal
関数が何もしないで終わってしまいます。
warn!
関数が警告文字を出力するためには、
出力用のhandler-bind
が必要となります。
出力用handler-bind
を用意する
これはプログラマーが明示的に用意するものではなく、
Common Lisp処理系があらかじめ用意しておくものです。
warning!
を出力するための、システムのコードを示します。
(handler-bind ((warning! (lambda (c) (format *error-output* "WARNING-TEST: ~A~%" c)))) ;; code )
上記の;; code
の部分で、eval-loop
であったり、
load
関数だったりが動作するわけです。
この要件から、warn!
関数は
(warn! "Hello")
としていたものは
(handler-bind ((warning! (lambda (c) (format *error-output* "WARNING-TEST: ~A~%" c)))) (warn! "Hello"))
みたいな感じになります。
warn!
関数がsignal
を起動するため、handler-bind
に制御が渡って、
format
関数によりWARNING-TEST
が出力されます。
一方、出力抑止である
(handler-bind ((warning! #'muffle-warning!)) (warn! "Hello"))
は、
(handler-bind ((warning! (lambda (c) (format *error-output* "WARNING-TEST: ~A~%" c)))) (handler-bind ((warning! #'muffle-warning!)) (warn! "Hello")))
となります。
warn!
関数はsignal
によりhandler-bind
を起動しますが、
内側のhandler-bind
に制御が渡り、mufffle-warning!
関数が呼び出されるため、
一番外側のhandler-bind
には制御が渡らず、WARNING-TEST
が出力されないのです。
この仕組みは、warning! condition
を途中で盗み見することができることを意味しています。
例えば次のコードを考えます。
(handler-bind ((warning! (lambda (c) (format *error-output* "<<<~A>>>~%" c)))) (warn! "Hello"))
出力結果
<<<Hello>>> WARNING-TEST: Hello
用意したhandler-bind
で一度condition
を受け取ってから、
format
で<<<Hello>>>
を出力させています。
condition
はそのまま外側に伝搬するので、
WARNING-TEST: Hello
も出力されるわけです。
それでは、これ以上伝搬させたくない場合はどうしたらいいでしょうか?
そんな時に使用するのがmuffle-warning
関数です。
(handler-bind ((warning! (lambda (c) (format *error-output* "<<<~A>>>~%" c) (muffle-warning! c)))) (warn! "Hello"))
出力結果
<<<Hello>>>
テストコード
説明で作成した自作のwarn!
関数を下記に示します。
(defun muffle-warning! (&optional condition) (let ((restart (find-restart 'muffle-warning! condition))) (if restart (invoke-restart restart) (error (make-condition 'control-error))))) (define-condition warning! (simple-condition) ()) (defmethod print-object ((instance warning!) stream) (let ((format (simple-condition-format-control instance)) (args (simple-condition-format-arguments instance))) (apply #'format stream format args))) (defun warn! (format &rest args) (restart-case (signal (make-condition 'warning! :format-control format :format-arguments args)) (muffle-warning! ()))) (defmacro progn-warning (&body body) `(handler-bind ((warning! (lambda (c) (format *error-output* "~&WARNING-TEST: ~A~%" c)))) ,@body)) ;; ;; main ;; (progn-warning (warn! "AAA") (warn! "BBB: ~A, ~A, ~A" 100 200 300) (handler-bind ((warning! #'muffle-warning!)) (warn! "CCC") (warn! "DDD")) (warn! "EEE"))
出力例
WARNING-TEST: AAA WARNING-TEST: BBB: 100, 200, 300 WARNING-TEST: EEE
clispだと2回出力されるのですが何故だろう。
無視します。
非正規化浮動小数の有効桁数を求める
C言語にてCommon Lispのfloat-precision
関数を作成したときのメモです。
非正規化の浮動小数とは、値が小さすぎるけど何とかして 浮動小数にて表現したいというギリギリの場合の方法です。
感覚としては、
1.23456 × 10の-103乗
みたいに表現したいけど指数-103は小さすぎて無理だから
0.00012 × 10の-99乗
のように表すのが非正規化です。
上記の例の場合、正規化の場合は1.23456
という値から
6桁まで表現できているのがわかりますが、
非正規化の場合は0.000
に押されて、12
という2桁にまで減ってしまっています。
この「2桁」というのをC言語の世界で求める方法を説明します。
一応注意書きしておくと、上記の例は基数が10になっていますが、 実際はだいたい2なので注意してください。
まず普通の時、つまり正規化された浮動小数の桁数は、
C言語のfloat.h
にdefineが用意されています。
FLT_MANT_DIG 24 DBL_MANT_DIG 53 LDBL_MANT_DIG 53, 64, 113とか?
もし浮動小数が正規化の場合は、何も考えずに上記の値を返却します。
では正規化か非正規化かを調べるにはどうしたらいいでしょうか。
思いつく方法は大小を比較することですが、
C99にはfpclassify
という関数が用意されています。
例えばdouble
を対象にした場合、
fpclassifyの戻り値が FP_NORMAL なら DBL_MANT_DIG
とします。 浮動小数点の値がゼロの場合はどうしましょうか? ここでは有効桁数は0ということにします。
fpclassifyの戻り値が FP_ZEOR なら 0
さて問題はfpclassify
の戻り値がFP_SUBNORMAL
の場合、つまり非正規化数の場合です。
返却値から先に示すと、
fpclassifyの戻り値が FP_SUBNORMAL なら DBL_MANT_DIG - (DBL_MIN_EXP - 1) + ilogb(value);
となります。
DBL_MANT_DIG
はすでに説明しました。
ilogb(value)
は、指数部を整数で返却する関数です。
問題はDBL_MIN_EXP
なのです。
DBL_MIN_EXP
は何なのでしょうか。
検索すると2通りの意味が出てきます。
- 基数の
e
乗で表現できる正規化数の最小値e
- 基数の
e-1
乗で表現できる正規化数の最小値e
どっちなんだ。
つまり最小値そのものなのか、そうじゃなく最小値に1
を足した数なのか?
現物を調査すれば早いので結果を言いますと、e-1
です。
IEEE754では、double
型の正規化数の指数部の最小値は-1022
になります。
つまり-1023
に相当する値(binary値だと0になる)が非正規化数となります。
正規化数の最小値-1022
より、
e-1 = -1022 e = -1021
となるので、
DBL_MIN_EXP = -1021
なのです。
では本当にC言語ではそうなのか?
コンパイラの問題だということもあり得ますので、
ちゃんと調べるためには規格書を見るしかありません。
ここであっているのかな。
- C - Approved standards
- The latest publically available version of the C11 standard is the document WG14 N1570, dated 2011-04-12.
Page 31 minimum negative integer such that FLT_RADIX raised to one less than that power is a normalized floating-point number, emin FLT_MIN_EXP DBL_MIN_EXP LDBL_MIN_EXP
英語わからないです。 でもone less thanとは、「1を引いた数」という意味のようです。 なるほど、なぜ1を引いているのかは知りませんが、規格がそうなっているんですね。
ということなので、
FP_SUBNORMAL なら DBL_MANT_DIG - (DBL_MIN_EXP - 1) + ilogb(value);
になります。
C言語にて、double
のコードを示します。
#include <float.h> #include <math.h> int float_precision_double(double v) { switch (fpclassify(v)) { case FP_NORMAL: return DBL_MANT_DIG; case FP_ZERO: return 0; case FP_SUBNORMAL: return DBL_MANT_DIG - (DBL_MIN_EXP - 1) + ilogb(v); default: fprintf(stderr, "error.\n"); exit(1); return 0; } }
実行の例を示します。
float_precision_double(123.456); →53 float_precision_double(0x01.0000p-1021); →53 float_precision_double(0x01.0000p-1022); →53 float_precision_double(0x01.0000p-1023); →52 float_precision_double(0x01.0000p-1024); →51 float_precision_double(0); →0
nptでは、ここで使っています。
最後に。
「非正規化」という名称はいろいろ種類があるように思えます。
デノーマル数という人もいれば、C言語では「SUBNORMAL」と表現されています。
Wikipediaの「非正規化数」が正しいのかもしれません。
rationalizeのアルゴリズムが全く分からなかった
Common Lispの関数rationalize
の話題です。
関数rationalize
とは、浮動小数を分数に変換する関数です。
例えば浮動小数の0.1
を分数に直すと次のようになります。
* (rationalize 0.1) 1/10
実は浮動小数の0.1
を1/10
に変換するって難しい事なのです。
0.25
を1/4
に変換するなら簡単です。
0.25
は基数2の浮動小数では自然に表せる数ですから。
0.1
は違います。
浮動小数の0.1
というのは、基数が2の浮動小数では
正確に表すことができずに仮数が循環してしまうものとして有名です。
例えば十進数で10
を3
で割ると0.3333...
になるのと同じです。
0.1
を単純に分数へ変換するにはrational
関数を使います。
* (rational 0.1) 13421773/134217728
すごい値ですね。
でもちゃんと0.1
なんです。
* (coerce 13421773/134217728 'float) 0.1
関数rational
の実装は簡単です。
仮数×2の指数乗
をそのまま分数形式にして通分するだけですから。
しかし関数rationalize
は違います。
どうも浮動小数0.1
と誤差の範囲から適切なものを探し出しているようなのです。
ではどうやって実装されているのでしょうか。
考えてもわかりませんでした。
仕方がないのでcmuclとsbclのソースを見てみました。
cmuclのソースを見ると、1618行目から下記のように始まっています。
;;; RATIONALIZE -- Public ;;; ;;; The algorithm here is the method described in CLISP. Bruno Haible has ;;; graciously given permission to use this algorithm. He says, "You can use ;;; it, if you present the following explanation of the algorithm."
ヒュー!! カッコいい!!
たぶんcmuclさんがrationalize
のアルゴリズムを使用したいってお願いして、
clispさんがいいよって言ったのかなって思います。
Bruno Haibleという名前には見覚えがあります。
なにせclispを起動するときに毎回出てくるのですから。
いつもお世話になっております。
$ clisp i i i i i i i ooooo o ooooooo ooooo ooooo I I I I I I I 8 8 8 8 8 o 8 8 I \ `+' / I 8 8 8 8 8 8 \ `-+-' / 8 8 8 ooooo 8oooo `-__|__-' 8 8 8 8 8 | 8 o 8 8 o 8 8 ------+------ ooooo 8oooooo ooo8ooo ooooo 8 Welcome to GNU CLISP 2.49 (2010-07-07) <http://clisp.cons.org/> Copyright (c) Bruno Haible, Michael Stoll 1992, 1993 Copyright (c) Bruno Haible, Marcus Daniels 1994-1997 Copyright (c) Bruno Haible, Pierpaolo Bernardi, Sam Steingold 1998 Copyright (c) Bruno Haible, Sam Steingold 1999-2000 Copyright (c) Sam Steingold, Bruno Haible 2001-2010 Type :h and hit Enter for context help. [1]>
こういうところで、Lisp仲間がつながっているのですね。
clispのソースはここで確認できます。
cmucl, sbcl, clispのソースには、コメントにアルゴリズムが載っています。
そのアルゴリズムを見ていたのですが、全く分からない。
sbclのソースには「An introduction to number theory」を見ろと書かれており、
さらにharvard.eduのリンク切れのサイトが紹介されている。
何となく思うんですが、
分数の割り算でつまずいている自分が
理解するのは非常に困難なんじゃないか。
理解できないなら既存のソースをたどって作るしかありません。 cmuclのソースはPublic Domainとの事なので、ありがたく参考にさせていただきます。
nptではcmuclのlispコードをもとにC言語で実装しました。 - https://github.com/nptcl/npt/blob/v0.1.6/src/real_decode.c#L746
他の関数はC言語に極力逆らわないように実装しているのですが、
関数rationalize
に関しては、C言語の皮をかぶったlispになっております。
しかしちゃんと動作しているようなので安心しました。
closureには何が保存されるのか
closureとは関数にデータを保存するための仕組みです。
Lisp大好きな人はlambda
と一緒に多用します。
例えばこんな使い方をします。
(let ((x 0)) (setq *call* (lambda () (incf x)))) (funcall *call*) -> 1 (funcall *call*) -> 2 (funcall *call*) -> 3 (funcall *call*) -> 4
上記の例では、lambda
により生成された無名関数が、
変数x
をclosureにて保存しています。
let
で宣言された変数x
のスコープはとっくに終わっているにもかかわらず、
funcall
で関数を呼び出すと、変数x
は立派に役目を果たしているというわけです。
では本題ですが、closureは何を対象にどんなものを保存するのでしょうか。
lexical変数だよ!と思った人は正解です。
問題はそれだけなのかということ。
典型的な回答としては下記の通り。
- lexical変数
flet
,labels
関数
flet
, labels
は局所関数を作成するためのものです。
変数ではなく関数です。
例えばこんな感じ。
(flet ((call () :hello)) (setq *call* (lambda () (call)))) (funcall *call*) -> :HELLO
面倒な話題として、setf
関数もclosureの対象となります。
説明すると長くなるので手短にしますが、
defsetf
とdefine-setf-expander
で定義されたマクロ形式の方ではなく、
flet
, labels
で定義した(setf 名前)
形式の関数がclosureの保存対象となります。
何が違うんだって思われるかもしれませんが、なんか結構違うんです。
そのうち詳しく説明します。
例をあげます。
(flet (((setf aaa) (value cons) (rplaca cons (* value value)) value)) (setq *call* (lambda (c v) (setf (aaa c) v)))) (let ((c (cons 10 20))) (funcall *call* c 999) c) → (998001 . 20)
ではdefun
宣言による関数はclosureに保存されないのでしょうか?
実はされません。
下記に例をあげます。
(defun aaa () :hello) (setq *call* (lambda () (aaa))) (funcall *call*) →:HELLO (defun aaa () :zzz) (funcall *call*) →:ZZZ
defun
による大域関数はspecial
変数みたいな扱いなんですね。
ここまではどのCommon Lispでも通用するでしょう。
ここからは処理系に依存する話になるかもしれません。
話は大きく変わり、tagbody
の実装をしたときの話です。
tagbody
のスコープを脱出したgo
について考えていました。
tagbody
は動的エクステントなので、スコープの外で呼び出されたgo
は無効となります。
問題は「無効」という意味が、放っておいてもいいオブジェクト
という意味では「ない」ということです。
無効になったtagは、呼び出された時点でエラーになります。
処理系は「エラーになる」という所まで面倒見てやらなければいけないのです。
例えば次の通り。
(tagbody (format t "Hello~%") again (format t "Tagbody~%") (setq *call* (lambda () (go again)))) (funcall *call*) →エラー、againはすでに閉じられている
どうすればいいのか。
go
に結びつけられたtagを専用のオブジェクトにするのです。
tagには戻り先のtagbody
に加えて、有効無効を表すフラグを持たせます。
tagはtagbody
とgo
で共有されます。
もしtagbody
がスコープを抜けたら、持っているすべてのtagに無効フラグを立たせます。
もしgo
が無効のtagを呼び出したらエラーにします。
では、そのtagオブジェクトは誰がいくつ保有するのか。
最初の考えだと、tagオブジェクトはtagbody
のコードが唯一保有すること考えていました。
つまり、(tagbody ...)
という表記が現れたら、
そのtagbody
がtagオブジェクトを保有するのだという考えです。
でもこれだと、再帰呼出か、スレッドによる同時実行により破綻します。
オブジェクトが一個しかないのですから。
解決方法は、tagオブジェクトはtagbody
を実行するたびに生成するしかなかったのです。
感覚としてはこんな感じ。
(let ((again (make-tagobject%))) (tagbody% (format t "Hello~%") again (format t "Tagbody~%") (setq *call* (lambda () (go% again))) (close-tagobject% again))
となると、当然tagオブジェクトもclosureに保存する必要が出てきます。
話が長くなりましたが、結局はtagbody
のtagもclosureに入れるんだよっていうことです。
全く同じ話がblock
/return-from
にも当てはまります。
しかしcatch
/throw
は対象外。
catch
はspecial
変数みたいな扱いなので必要ありません。
結論は、closureには下記のオブジェクトが格納される。
- lexical変数
flet
,labels
関数flet
,labels
のsetf
関数tagbody
のtagblock
のsymbol
分数の割り算で余りを求める
分数の割り算で、余りを求める方法がわからなかったという内容です。
どうも小学校の算数の内容らしいですが、わからなかったのだから仕方がない。
あと「帯分数」というモノを初めて知りました。
習った記憶が全くない。
Common Lispでは、truncate
関数を使うことによって求められます。
では計算方法はどうなるのでしょうか。
例えば「100/3 ÷ 7/3」の商と余りはいくつになるでしょう。
素直に割り算すると100/7になるので、
「14 余り 2/7」が答えになりそうですが違います。
これだと割り算の結果を、ただ整数と分数の部分に分けているだけです。
つまりはこんな感じ。
* (truncate (/ 100/3 7/3)) 14 ; 2/7
正しい回答を出力してみます。
* (truncate 100/3 7/3) 14 ; 2/3
商の14までは合っていますが、余りが違います。 この余りをちゃんと求めるには、商と余りの関係を式で表す必要があります。
割り算
a ÷ b = 商...余り
のとき、書き換えると
余り = a - b×商
なので、
余り = 100/3 - (7/3)×14
計算をすると余り = 2/3
となります。
この分数の余りの求め方はfloor
, ceiling
, truncate
, round
,
mod
, rem
のすべてで使っています。
当たり前と言えばそうなんですが、実装面だと割り算で分数を扱うときに
初めて上記の計算が必要になります。
分数以外の型であるfixnum
, bignum
, float
は別の方法で求めるので、
この計算を自分でしなくてもいいというのがなんとなく不思議な感じがします。
Common LispとC言語で虚数部のゼロ符号に違いが出る
ずっと悩んでいたことがわかったのでメモします。
例えば下記の関数
log(1-z)
のzに10を入れると、Common LispとC言語では結果が違います。
(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
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)
が含まれるので、虚数部ゼロの符号に差異が生じるというわけです。
何が正しいのか全くわからん。
C言語にて関数csinlが使えない場合がある
数学のsin関数に複素数って渡せるんですね。 全然知りませんでした。
C言語のc99標準ライブラリにも用意されているようです。
具体的には下記の関数。
csinf, csin, csinl, sin(tgmath)
問題は、clangのlong double complexを引数とする関数がいくつか使えない事です。
例えばclangではcsinlを使おうとするとエラーになります。
clangというとFreeBSDの標準コンパイラですが、
manコマンドのマニュアルにも載ってなようなので、
少なくともFreeBSD開発側は認識できているようです。
確認した環境は下記の通り。
% freebsd-version 12.0-RELEASE-p3 % clang --version FreeBSD clang version 6.0.1 (tags/RELEASE_601/final 335540) (based on LLVM 6.0.1) Target: x86_64-unknown-freebsd12.0 Thread model: posix InstalledDir: /usr/bin
モジュール自体が存在しないので、csinlがあるとコンパイルが通りません。 tgmath.hを利用してsin関数にlong double complexを渡しても、 あれは結局マクロでcsinlにするだけなのでエラーになります。
そうなんですか、無理ですか。
完全に他人任せで実装されるまで大人しく待ってようと思います。
しかし妥協すれば解決策はいくつかあると思います。
私の大のお気に入りは、下記の極めて最低な方法。
#ifdef __clang__ #define csinl csin #endif
#include <stdio.h> #include <cmath> #include <complex> int main() { std::complex<float> f(1.2f, 3.4f); std::complex<double> d(1.2, 3.4); std::complex<long double> l(1.2L, 3.4L); f = std::sin(f); d = std::sin(d); l = std::sin(l); printf("%30.25e, %30.25e\n", real(f), imag(f)); printf("%30.25e, %30.25e\n", real(d), imag(d)); printf("%30.25Le, %30.25Le\n", real(l), imag(l)); return 0; }
実行結果
$ clang++ test.cpp $ ./a.out 1.3979410171508789062500000e+01, 5.4228153228759765625000000e+00 1.3979408806017994848502894e+01, 5.4228154724634016758955113e+00 1.3979408806017997938912767e+01, 5.4228154724634012448167275e+00
うまく行っているように見えます。
なぜC++でうまく行ってC99じゃダメなのか。
他には、clangではなくgccを利用する方法。
2月くらいのGentoo Linuxで確認しましたが、こちらは問題ありません。
そしてあまりお勧めできないのは自分でcsinlを実装する方法です。
正直csinl程度なら何とかなりそうなもんですけど、
arcsinやらになると、数値計算に自信がないならやめた方がいいかもしれません。
もしかして自分が勘違いしているだけでcsinlはc99標準じゃない? 関数clogなんて丸ごと存在しないんですよね。 わざわざ/usr/includeを確認してしまった。