libm
Windows10に早くしろしろウィルスに感染したWindows7機、約5時間待ってみたけど、 さっっぱり更新が始まらない。観念して個別にインストールしようと思って、 パッチをDLし始めた。
DLスピードはめちゃ遅い。60kbpsぐらいしか出ない。サーバー側は未だにADSL端子に なってるんではなかろうか？ それとも７ユーザーに嫌がらせの意味を込めて、 優先度を最低ランクの端子にしてないか？
そんなこんなでDL中に、ふとアップデートの画面を見たら、やっとDLできまっせ案内が 来てた。開始してから5時間半ぶりですよ。そして、総DL容量は86Mですって。ここからも 忍耐の時間が始まるのね。もー、やだやだ。付き合っていられないな。
libm
前回、数学で最も美しい式と言われる、オイラーさんが見つけた式を、各種言語で やってみた。その結果、オイラーさんが（きっと）現代コンピュータ言語を嘆いて いるに違いないと思った。
虚部の項はZEROになるはずなのに、若干の値が残っている。原因はきっとlibm.so に有るだろうと推測。eのべき乗を求める関数 exp は、超越関数。実行スピードを 稼ぐ為に、あの糞石メーカー製作の i387が使われているに違いない。
ならば、どんなOSを持ってこようと、x86で動くOS上なら、同じ結果になるはず。 暇なんで、検証してみる。
検証方法は、gaucheを基本として、定数の pi は、 2 * acos(0) で、計算させる 事にした。これなら、定数の定義桁数によらない値を使ってくれると思ったから。
$ gosh gosh> (exp (* 0+1i 2 (acos 0))) -1.0+1.2246063538223773e-16i
実部の表示は、1 だったり、1. だったりするけど、見やすいように、1.0 とした。
FreeBSD -1.0+1.2246063538223773e-16i ; gauche OpenBSD -1.0+1.2246063538223773e-16i ; gambit NetBSD -1.0+1.2246063538223773e-16j ; python Fedora -1.0+1.2246467991473532e-16i ; gauche Ubuntsu -1.0+1.2246467991473532e-16i ; gauche iOS9 -1.0+1.2246467991473532e-16i ; gambit iOS9 -1.0+1.22464679915e-16j ; pythonista
何やら、BSD vs. Linux と言う構図だな。これはもうソース嫁の世界に突入必至か。 まてまて、オイラーは、ARMな石も持ってるぞ。この石なら、どう転んでも、i387は 使えまえ。iOS9にも参入して貰った。もう、ipadを使い始めて3年経過したけど、 こんな事したのは初めてだ。
iOS族もLinuxに右に倣えしてるな。アプルはlibmみたいのを共通基盤として提供してる？ それはちょっと考えづらいな。libmなんてわざわざ指定しないとリンク出来ない人気の なさ。と言うか特殊なもの。
あれ、画像関係でsinとかcosが多用されるけど、どうなってるの？ そんなものは、 画像データ生成・加工機のGPUに丸投げさ。この辺を90度回転してって下請けに出す だけ。だから、libmは不要と思われ。
ユーザーが使っていいのは、画面の描画とか通信とかの派手系ぐらいだろう。 用意したとしても、ゲームで多用されるであろう、乱数の発生ぐらいか。ああ、乱数は セキュリティの要になるから、絶対アプル製だろうね。
以前、appの構成をちら見したけど、appが使うライブラリィ類は、 appの中に閉じていたように思う。って事は、app作者が自前で用意するしかない。
人気はLinux系って事で、GNUのソースをコンパイルしてappに直付けしたか、ライブラリィ で持つ事にしたのだろう。だから、指紋が残ったとな。
gambitのソースは公開済み。pythonistaって言うpython系は、ソースが公開されてるのかな？ 若し未公開なら、GPLに汚染してる証拠を提示して、公開を迫ってみようか。
例によってLinux系のソースは入れてないので、OpenBSDのlibmを紐解いてみますかな。
複素数演算の探検用ソース
libmを探検する前に、探検用のツールを作っておく。goshのコードをC語に置き換えた やつを書く。
#include <stdio.h>
#include <complex.h>
#include <math.h>
int main(void){
long double _Complex z, x;
x = I * 2 * acosl(0);
z = cexpl(x);
printf("real=%g\timag=%38.36g\n", creal(z),cimag(z));
return 0;
}
complex.hを使うと、虚数単位(0+1i)を大文字のIとして使えるようになる。また、πは、M_PIとして 定義されてるので、それを利用してもいいのだけど、goshに習って自前で導出した。 全体は、long doubleで無駄に精度よく計算してみた。複素数の関数の説明は、 標準ライブラリ(新規)が詳しい。C99規格 なんですねぇ。C語も地味に進化してるな。
違うって、複素数と言えば昔からFortran。この牙城を崩すには、C陣営も複素数について 早く規格化する事が求められた。それで、21世紀になる前に、どうしても決める必要が 有ったのさ。だから、誇らしげにC99って言ってる。願わくば、21世紀はC語が覇者と たらん事を。ああ、それなのにそれなのに、Javascriptがしゃしゃり出てきおって。 Cフラフラはしぶとく、改定を重ねているな。
使い方は、
[ob: z]$ cc -g t.c -lm [ob: z]$ ./a.out real=-1 imag=-5.01655761266833226942395050995492469e-20
演算精度を高めたら、誤差項が小さくなったぞ。出来上がった a.outをobjdump -S して、コードを覗いてみる。
z = cexpl(x); ba3: 8d 75 88 lea 0xffffff88(%ebp),%esi ba6: 8b 45 a8 mov 0xffffffa8(%ebp),%eax ba9: 8b 55 ac mov 0xffffffac(%ebp),%edx bac: 8b 4d b0 mov 0xffffffb0(%ebp),%ecx baf: 89 44 24 04 mov %eax,0x4(%esp) bb3: 89 54 24 08 mov %edx,0x8(%esp) bb7: 89 4c 24 0c mov %ecx,0xc(%esp) bbb: 8b 45 b4 mov 0xffffffb4(%ebp),%eax bbe: 8b 55 b8 mov 0xffffffb8(%ebp),%edx bc1: 8b 4d bc mov 0xffffffbc(%ebp),%ecx bc4: 89 44 24 10 mov %eax,0x10(%esp) bc8: 89 54 24 14 mov %edx,0x14(%esp) bcc: 89 4c 24 18 mov %ecx,0x18(%esp) bd0: 89 34 24 mov %esi,(%esp) bd3: e8 f8 fb ff ff call 7d0 <__init+0x50> bd8: 83 ec 04 sub $0x4,%esp bdb: 8b 45 88 mov 0xffffff88(%ebp),%eax bde: 8b 55 8c mov 0xffffff8c(%ebp),%edx be1: 8b 4d 90 mov 0xffffff90(%ebp),%ecx be4: 89 45 c8 mov %eax,0xffffffc8(%ebp) be7: 89 55 cc mov %edx,0xffffffcc(%ebp) bea: 89 4d d0 mov %ecx,0xffffffd0(%ebp) bed: 8b 45 94 mov 0xffffff94(%ebp),%eax bf0: 8b 55 98 mov 0xffffff98(%ebp),%edx bf3: 8b 4d 9c mov 0xffffff9c(%ebp),%ecx bf6: 89 45 d4 mov %eax,0xffffffd4(%ebp) bf9: 89 55 d8 mov %edx,0xffffffd8(%ebp) bfc: 89 4d dc mov %ecx,0xffffffdc(%ebp)
盛大にmov命令が並んでいるな。longのdoubleを32Bitマシンでやると、こういう塩梅に なるのか。早く来い来い64Bitマシン。0xbd3が、関数cexplへの入り口のようだ。
ちょいとgdbで追ってみる。
Breakpoint 1, cexpl (z=0 + 3.1415926535897932385128089594061862 * I)
at /usr/src/lib/libm/src/s_cexpl.c:66
66 r = expl(creall(z));
:
main () at t.c:11
11 printf("real=%g\timag=%38.36g\n", creal(z),cimag(z));
もう、mainに戻ってきた。滞空時間はわずか。でも、ここ見ろって教えられたから 良しとしよう。
60long double complex
61cexpl(long double complex z)
62{
63 long double complex w;
64 long double r;
65
66 r = expl(creall(z));
67 w = r * cosl(cimagl(z)) + (r * sinl(cimagl(z))) * I;
68 return (w);
69}
もう全く数学の世界です。複素数が実世界にひょいと顔を出す。で、このdir近辺を見ると、
[ob: libm]$ ls Makefile man/ obj@ src/ arch/ noieee_src/ shlib_version [ob: libm]$ ls arch/i387/ e_acos.S e_sqrtf.S s_floor.S s_rint.S e_asin.S e_sqrtl.S s_floorf.S s_rintf.S e_atan2.S fenv.c s_ilogb.S s_scalbnf.S e_atan2f.S invtrig.c s_ilogbf.S s_significand.S e_exp.S s_atan.S s_llrint.S s_significandf.S e_fmod.S s_atanf.S s_llrintf.S s_sin.S e_log.S s_ceil.S s_log1p.S s_sinf.S e_log10.S s_ceilf.S s_log1pf.S s_tan.S e_remainder.S s_copysign.S s_logb.S s_tanf.S e_remainderf.S s_copysignf.S s_logbf.S e_scalb.S s_cos.S s_lrint.S e_sqrt.S s_cosf.S s_lrintf.S
一応、i387用のコードも置いてある。ちょいと、i387/e_exp.Sを見ておくか。
/* e^x = 2^(x * log2(e)) */
ENTRY(exp)
:
fldl2e
fmulp /* x * log2(e) */
fst %st(1)
frndint /* int(x * log2(e)) */
fst %st(2)
fsubrp /* fract(x * log2(e)) */
f2xm1 /* 2^(fract(x * log2(e))) - 1 */
fld1
faddp /* 2^(fract(x * log2(e))) */
fscale /* e^x */
:
糞石で直接expは計算出来ないから、別の方法を取るとな。これが、i387の実体。 昔はスピードを稼ぐために持て囃されたけど、アメリカの団体 IEEE が、押す ルーチンを使えやゴリャと、押し付けが有った模様。
そのせいで、i387も、反IEEEなルーチンも、過去の遺物として、ひっそりと置いて あるって事だな。米国のごり押し、恐るべしが結論ですかね。noieee_srcの中を みると、
* Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved.
こんなクレジットが入っている。学が産業界に押されたって事か。
でも、普通に考えると、不要なソースを入れて（保持）おくはずが無い。どういう風に コンパイルされてるか調べれば良いだろう。Makefileを読むのも手だけど、擬似実行して ログを眺めるのが楽かな。
[ob: libm]$ make -n | tee ~/MLOG
echo "cc -g -I/usr/src/lib/libm/src -I/usr/src/lib/libm/src/ld80 -c /usr/src/lib/libm/arch/i387/e_acos.S -o e_acos.o" cc -g -I/usr/src/lib/libm/src -I/usr/src/lib/libm/src/ld80 -c /usr/src/lib/libm/arch/i387/e_acos.S -o e_acos.o.o ld -X -r e_acos.o.o -o e_acos.o rm -f e_acos.o.o : echo building standard m library rm -f libm.a ar cq libm.a `lorder b_exp__D.o b_log__D.o b_tgamma.o e_acos.o e_acosf.o .... : echo "cc -g -I/usr/src/lib/libm/src -I/usr/src/lib/libm/src/ld80 -c -DPROF /usr/src/lib/libm/arch/i387/e_remainder.S -o e_remainder.po" cc -g -I/usr/src/lib/libm/src -I/usr/src/lib/libm/src/ld80 -c -DPROF /usr/src/lib/libm/arch/i387/e_remainder.S -o e_remainder.po.o ld -X -r e_remainder.po.o -o e_remainder.po rm -f e_remainder.po.o : echo building profiled m library rm -f libm_p.a ar cq libm_p.a `lorder b_exp__D.po b_log__D.po b_tgamma.po e_acos.po ... : echo "cc -g -I/usr/src/lib/libm/src -I/usr/src/lib/libm/src/ld80 -c -fpic /usr/src/lib/libm/arch/i387/e_atan2f.S -o e_atan2f.so" cc -g -I/usr/src/lib/libm/src -I/usr/src/lib/libm/src/ld80 -c -fpic /usr/src/lib/libm/arch/i387/e_atan2f.S -o e_atan2f.so.o ld -X -r e_atan2f.so.o -o e_atan2f.so rm -f e_atan2f.so.o : echo building shared m library \(version 9.0\) rm -f libm.so.9.0 cc -shared -fpic -o libm.so.9.0 `lorder b_exp__D.so b_log__D.so ...
三回もコンパイルしてる。静的用、それのプロファイル版、そしてsoファイルとな。 ld80ってのをincludeのdirに指定してi387の下のやつをコンパイルしてるけど、これって、 出来たらインテル特製の80Bit浮動点機構を使ってねって事かしら。id128ってのも有るけど、 これは64Bitに進化した石用かな。
i387(FPU) の現場
上で作った簡単アプリは、精度が128Bit用。これじゃ、糞石のハードによる演算機構は 及びじゃない。64Bit精度にしたら、i387の出番も有るだろう、と言う目論見で、 やってみる。
出来たa.outを逆アセしてみると、
z = cexp(x); b57: 8d 45 b8 lea 0xffffffb8(%ebp),%eax b5a: dd 45 d0 fldl 0xffffffd0(%ebp) b5d: dd 5c 24 04 fstpl 0x4(%esp) b61: dd 45 d8 fldl 0xffffffd8(%ebp) b64: dd 5c 24 0c fstpl 0xc(%esp) b68: 89 04 24 mov %eax,(%esp) b6b: e8 40 fc ff ff call 7b0 <__init+0x30> b70: 83 ec 04 sub $0x4,%esp b73: dd 45 b8 fldl 0xffffffb8(%ebp) b76: dd 5d e0 fstpl 0xffffffe0(%ebp) b79: dd 45 c0 fldl 0xffffffc0(%ebp) b7c: dd 5d e8 fstpl 0xffffffe8(%ebp)
今までのmov一辺倒が、劇的に変わった。ニューモニックの頭に付いている、f は、 i387専用の命令だろう。
i387な石にBPを置けるか知らないけど、やってみる。
(gdb) b /usr/src/lib/libm/arch/i387/e_exp.S:52
No source file named /usr/src/lib/libm/arch/i387/e_exp.S.
Make breakpoint pending on future shared library load? (y or [n]) y
Breakpoint 1 (/usr/src/lib/libm/arch/i387/e_exp.S:52) pending.
(gdb) run
Starting program: /home/sakae/z/a.out
Breakpoint 1, exp () at /usr/src/lib/libm/arch/i387/e_exp.S:52
52 fldl 4(%esp)
(gdb) bt
#0 exp () at /usr/src/lib/libm/arch/i387/e_exp.S:52
#1 0x02b19f4b in cexp (z=0 + 3.1415926535897931 * I)
at /usr/src/lib/libm/src/s_cexp.c:68
#2 0x14d12b70 in main () at d.c:9
:
(gdb) si
77 fstp %st(1)
(gdb) info float
R7: Empty 0x3fffb8aa3b295c17f0bc
R6: Empty 0x80000000000000000000
R5: Empty 0x80000000000000000000
R4: Empty 0x00000000000000000000
R3: Empty 0x00000000000000000000
R2: Zero 0x00000000000000000000 +0
=>R1: Valid 0x3fff8000000000000000 +1
R0: Empty 0x3fff8000000000000000
Status Word: 0x0900 C0
TOP: 1
Control Word: 0x037f IM DM ZM OM UM PM
PC: Extended Precision (64-bits)
RC: Round to nearest
Tag Word: 0xffd3
Instruction Pointer: 0x00:0xfc21b6d4
Operand Pointer: 0x00:0x00000000
Opcode: 0xd9fd
へぇー、レジスターが８本と、若干の付属レジスターが有るとな。64bitの精度も 許しているのか。続けて、sinも追ってみる。
69 w = r * cos (y) + r * sin (y) * I;
(gdb) b /usr/src/lib/libm/arch/i387/s_sin.S:10
Breakpoint 2 at 0x2b31104: file /usr/src/lib/libm/arch/i387/s_sin.S, line 10.
(gdb) c
Continuing.
Breakpoint 2, sin () at /usr/src/lib/libm/arch/i387/s_sin.S:10
10 fldl 4(%esp)
(gdb) bt
#0 sin () at /usr/src/lib/libm/arch/i387/s_sin.S:10
#1 0x0ea88f69 in cexp (z=0 + 3.1415926535897931 * I)
at /usr/src/lib/libm/src/s_cexp.c:69
#2 0x168e3b70 in main () at d.c:9
(gdb) si
11 fsin
(gdb) si
12 fnstsw %ax
(gdb) si
13 andw $0x400,%ax
(gdb) info float
R7: Empty 0x3fffb8aa3b295c17f0bc
R6: Empty 0x80000000000000000000
R5: Empty 0x80000000000000000000
R4: Empty 0x00000000000000000000
R3: Empty 0x00000000000000000000
=>R2: Valid 0x3fca8d30000000000000 +1.224606353822377258e-16
R1: Empty 0x3fff8000000000000000
R0: Empty 0x3fff8000000000000000
Status Word: 0x1320 PE C0 C1
TOP: 2
Control Word: 0x037f IM DM ZM OM UM PM
PC: Extended Precision (64-bits)
RC: Round to nearest
Tag Word: 0xffcf
Instruction Pointer: 0x00:0xfc21b6d8
Operand Pointer: 0x00:0x00000000
Opcode: 0xd9fe
(gdb) c
Continuing.
real=-1 imag=1.22460635382237725821141793858259916e-16
[Inferior 1 (process 22731) exited normally]
/usr/src/lib/libm/src/s_cexp.c の sin 呼び出しが、いつの間にか、 /usr/src/lib/libm/arch/i387/s_sin.Sのマシン語呼び出しに展開されている。 その裏に、OS側が、色々やってるんだけどね。gdbで追っていくと、いやになる程、 色々なルーチンを走り回っている。この機構って、ダイナミックリンクの、ジャンプ テーブル書き換えの技が使われているのだろうね。
Linux系はどうよ？
比べてみる。
[sakae@fedora z]$ ./a.out real=-1 imag=-5.01655761266833226942395050995492469e-20
sakae@arm:~/z$ uname -a Linux arm 3.16.0-4-versatile #1 Debian 3.16.7-ckt9-3~deb8u1 (2015-04-24) armv5tejl GNU/Linux sakae@arm:~/z$ ./a.out real=-1 imag=1.2246467991473532071737640294583966e-16
armなやつは、long doubleをサポートしていないな。昔な石なんで、ソフトのサポートも サボったんかな。一応、debuggerに通してみた。
(gdb) bt
#0 __sincos (x=3.1415926535897931, sinx=0xbefff5f0, sinx@entry=0xbefff5d8,
cosx=0xbefff5f8, cosx@entry=0xbefff5e0)
at ../sysdeps/ieee754/dbl-64/s_sincos.c:35
#1 0xb6f67000 in __cexp (x=0 + 3.1415926535897931 * I) at s_cexp.c:44
#2 0x0001054c in main () at t.c:9
debian系はソースを取ってこないと参照出来ない。fedoraで無理して参照しようと すると、glibc-debuginfoを gdbのお供にいかがでしょうかって勧められて、それを入れると一応gdb上から 便利にソースが見られる。で、libmって、しっかり GNU Cの仲間って書いてあったぞ。
/* Return value of complex exponential function for double complex value.
Copyright (C) 1997-2015 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
それはいいんだけど、いきなりコードが始まっていて、数学的な根拠が何も 書いてない。なんだかなあ。
BSDはどうよ？
対抗上、今度はFreeBSDのソースを三田。どれでもいいんだけど、/usr/src/lib/msun/src/k_cos.c と言う、 cosの計算の基本となるものにした。dirツリーがlibmじゃなくて、msunになってるのは、 vimの作者であるビル・ジョイ が居た会社、Sunに敬意を表しているのかな。
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
:
* __kernel_cos( x, y )
* kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
*
* Algorithm
* 1. Since cos(-x) = cos(x), we need only to consider positive x.
* 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
* 3. cos(x) is approximated by a polynomial of degree 14 on
* [0,pi/4]
* 4 14
* cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
* where the remez error is
*
* | 2 4 6 8 10 12 14 | -58
* |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
* | |
*
* 4 6 8 10 12 14
* 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
* cos(x) ~ 1 - x*x/2 + r
* since cos(x+y) ~ cos(x) - sin(x)*y
* ~ cos(x) - x*y,
* a correction term is necessary in cos(x) and hence
* cos(x+y) = 1 - (x*x/2 - (r - x*y))
* For better accuracy, rearrange to
* cos(x+y) ~ w + (tmp + (r-x*y))
* where w = 1 - x*x/2 and tmp is a tiny correction term
* (1 - x*x/2 == w + tmp exactly in infinite precision).
:
whereとかletなんて単語が出てきて、何となくHaskellっぽい雰囲気で説明が 続いている。ちゃんと根拠が有るのは、産業界では必須な事ですよ。
ここまでは良いんだけど、NetBSDで試してみたら、cexplは受け付けてくれるものの 精度はdoubleでの計算しか出来ずのズルモードが発動しました。 FreeBSDでは、コンパイルの段階で、complex.hをincludeせいや、あるいはテメーで 定義しろやと抜かしおって、バイナリーを作ってくれず。同じ血を分けたBSD３兄弟と 言えども、性格が出てて面白い。
Intel 8087
FreeBSD上の呼称は i387 一般的にはFPU。もっと一般化すると、コ・プロセッサーに なるかな。armとかmipsな石だと、コプロとか言うね。
Intel 8087(日本語版の概略説明)
Intel 8087 (英語版、詳しい説明有り)