fdlibm
Table of Contents
maxima
Maxima is a Common Lisp implementation of MIT's Macsyma system for computer based algebra.
sbcl/gcl/ecl とかのコモン・リスプ上に構築された数式処理システム。これ で何をやろうとしてる？ 前回やったawkの組み込み関数である、三角関数を探 てみたいんだ。
なにが使えるかと言うと、テーラー展開。使用方法は、
taylor(f(x), x, a, n) -> f(x) を x=a の周りで n 次の項までテイラー展開する
こんな具合だ。下記の例だと、sin関数は、別の表現ができますよって事。 atanの方は、分母が素直だな。sinの方は、何となく階乗かな。どんな数列か 当てるには、それなりの経験が必要っぽい。
(%i1) taylor(sin(x), x, 0, 10);
3 5 7 9
x x x x
(%o1)/T/ x - -- + --- - ---- + ------ + . . .
6 120 5040 362880
(%i2) taylor(atan(x), x, 0, 10);
3 5 7 9
x x x x
(%o2)/T/ x - -- + -- - -- + -- + . . .
3 5 7 9
x=a の周りで って、どゆ事？ ひょっとして接線とかと関係ある？
やっぱり、思った通りだ。高次微分ができる関数ならば多項式で近似できる。
こちらは、格調高い説明。なんとなく、なんとなくの感がするな。
FPUと言えども、sinの算出方法を知っている訳ではない。せいぜい計算できる のは、加算と乗算ぐらいだ。それを使ってsinとかを計算しようとしたら、積 和の世界まで分解する必要がある。それで数学の一分野、テーラー展開のおで ましだ。
数学って、発見なの？ それとも発明なの？ その両方があって発展してきたん だろうな。
softfpu
まずは、整数じゃない世界。数学的に言うと実数をコンピューターにも持ち込 む方法。
Berkeley SoftFloatを使って浮動小数点数演算をソフトウェアエミュレートで置き換える(x8664)
Rustのsoftfloatライブラリ比較 ここで紹介されてる rug っての、任意精度 で使えるらしい。と言うことは、AI用なんで精度は必要ありません。乗算と加 算ができれば、行列には十分ですから、とかなるな。
qemu
ひょっとして、こいつに参考になりそうなコードが有るかな？
qemu/target/i386/tcg/fpu_helper.c
void helper_fcos(CPUX86State *env)
{
double fptemp = floatx80_to_double(env, ST0);
if ((fptemp > MAXTAN) || (fptemp < -MAXTAN)) {
env->fpus |= 0x400;
} else {
ST0 = double_to_floatx80(env, cos(fptemp));
env->fpus &= ~0x400; /* C2 <-- 0 */
/* the above code is for |arg| < 2**63 only */
}
}
たとえば、こんなのが見付かるけど、これって、土台はインテルのi387に依存 してないか。もっとクリーン・コンピュータが欲しいぞ。
gdb/sim
gdbにもシュミレータが付属してたよね。探してみると、楽しいのが出てきた。
gdb-13.2/sim/common/sim-fpu.c
/* @(#)e_sqrt.c 5.1 93/09/24 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* __ieee754_sqrt(x)
ただ、sqrtだけなんだよな。
fdlibm
と言う事で、前回発見した、サン マイクロシステムからの遺産です。
https://netlib.org/fdlibm/ ieee754と言うベーシックな規格やら、派生？の 規格も含めて実装されている。
ダウンロードは、 wget -r https://netlib.org/fdlibm/ したんだけど、curlの場合はどうすればいいの？
簡単なテストコードを作成。
#include "fdlibm.h"
#include <stdio.h>
int main(){
printf("%.10g\n" , 4 * atan( 1.0 ) );
}
普通なら、これでOK. ああ、その前に gmake するのをお忘れなく。
[sakae@deb fdlibm]$ cc zmain.c libm.a [sakae@deb fdlibm]$ ./a.out 3.141592654
コンパイルした時の残骸が残っていたなら、これでもOK.
vbox$ cc zmain.c s_atan.o vbox$ ./a.out 3.141592654
gdbが使える様にコンパイルしとくと
vbox$ gdb -q a.out Reading symbols from a.out... (gdb) b atan Breakpoint 1 at 0x1860: file s_atan.c, line 95. (gdb) r Starting program: /tmp/fdlibm/a.out Breakpoint 1, atan (x=1) at s_atan.c:95 95 hx = __HI(x); (gdb) bt #0 atan (x=1) at s_atan.c:95 #1 0x1a06780e in main () at zmain.c:5
以下、ステップ実行での主要な流れ。
124 z = x*x; 125 w = z*z; 127 s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10]))))); 128 s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9])))); 129 if (id<0) return x - x*(s1+s2); 131 z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); 132 return (hx<0)? -z:z; 134 } (gdb) 3.141592654 main () at zmain.c:6 6 }
at debian
同じ事をdebianでやってみると、サン・マイクロの遺産と合体しなくても、こ の通りに、実行できちゃった。
[sakae@deb fdlibm]$ gcc zmain.c
[sakae@deb fdlibm]$ ./a.out
3.141592654
[sakae@deb fdlibm]$ ldd a.out
linux-gate.so.1 (0xb7f9a000)
libc.so.6 => /lib/i386-linux-gnu/libc.so.6 (0xb7d8a000)
/lib/ld-linux.so.2 (0xb7f9c000)
ひょっとして、お前、余計な事してないか？ たとえば、libc.so.6に該当する コードが入れてあるとか。
前回、debianでgdb/awkしてた時、ソースが無いなんて言ってたけど、それって glibcの一員だったりしないか。忖度は国府だけで十分(それは、いけない事だ けど)と思うぞ。オイラーはクリーンな環境が欲しいんだ。忖度まみれは、御 免被る。
まて、そんなの、いいがかりだろう。時代錯誤もはなはだしいぞ。今やFPUは CPUの中に吸収されちゃって分離不能。だったら、わざわざFPUを叩くコード (libmの事ね)を分離する必要なんてない。だったらlibcに吸収しちゃって、何 も考えずに使える方が便利。こういうご時世なんだな。
組み込みシステムとかFPUの無い世界で開発をする場合(ごく少数の人だろうけ ど)は、このglibcを自動リンクしないオプションがgccに用意されてるんだな。
at FreeBSD
FreeBSDには色々な版のgccが導入されているので、確認してみる。まずは、生 粋のclangで、正しい挙動の確認から。
[sakae@fb /tmp/fdlibm]$ cc zmain.c ld: error: undefined symbol: atan >>> referenced by zmain.c >>> /tmp/zmain-bf996f.o:(main) cc: error: linker command failed with exit code 1 (use -v to see invocation) [sakae@fb /tmp/fdlibm]$ gcc12 zmain.c [sakae@fb /tmp/fdlibm]$ ./a.out 3.141592654
gcc10,gcc11もdebian風に動いた。
maximaで追跡
判りやすい様に、0.5 = sin(x) を想定しよう。さて、このxは何度になるでしょ うか？ 何度じゃなくて難度と言うのは無しよ。そんなの暗算しょ。
正三角形を思いうかべる。どっしりとしたやつね。その頂点から垂線を下す。 直角三角形が2つできた。おのおのの三角形の頂点側の角度は30度になる。 それに対向する辺(底辺)の長さは、正三角形の 1/2 だ。よって、sinの定義から、辺 /斜辺は、0.5。その時の角度は30度。ラジアンで言うと pi/6となる。
(%i1) taylor(sin(x), x, 0, 7);
3 5 7
x x x
(%o1)/T/ x - -- + --- - ---- + . . .
6 120 5040
(%i2) define(ms7(x), %o1);
3 5 7
x x x
(%o2)/T/ ms7(x) := x - -- + --- - ---- + . . .
6 120 5040
(%i3) ms7(%pi / 6);
7 5 3
%pi - 1512 %pi + 1088640 %pi - 235146240 %pi
(%o3)/R/ - -----------------------------------------------
1410877440
(%i4) float(%o3);
(%o4) 0.49999999186902316 ; n = 7
(%o8) 0.5000000000202799 ; n = 9
(%o12) 0.4999999999999643 ; n = 11
(%o16) 0.4999999999999999 ; n = 13
(%i17) fpprec:40;
(%o17) 40
(%i20) bfloat(%o20);
(%o20) 5.000000000000000465631247757710526025216b-1
(%i1)で、7項までテーラー展開。(%i2)で、その結果の%o1 を右辺とする関数 ms7を定義。(%i3)で、計算。結果は、定数の%piを含めた式になる。だって、 展開しちゃうと、正確な表現が消失しますからね。
でも、現実にどれぐらいの値になるか知りたいはず。そんな時は、数値に変換 してねって指示すればよい。そんな要領で展開次数を増加させてみた。
なお、無限精度の計算をしたい場合、何桁までって指示と、それ用の計算指示 をすればよい。
で、この結果から、項数を増やすと、より正解に近付く。あれと一緒だな。 そう、sin波を複数合成すると、方形波でも鋸歯状波でも自在に作成できる。 より正解を目指すなた、沢山の波が必要ってね。
そう考えると、この式こそ、生成多項式を名乗って欲しいぞ。sin値が、それ とは全く異なる方法で計算されるんですから。QRコードの誤り検出で、生成多 項式が出てきたけど、オイラーには意味不な語句であった。
sinの追跡
本物の math.h には、こんな定義がされているので、利用させてもらう。
#define M_PI ((double)3.14159265358979323846) /* pi */ #define M_PI_2 ((double)1.57079632679489661923) /* pi/2 */ #define M_PI_4 ((double)0.78539816339744830962) /* pi/4 */
sin(45度)の計算
#include "fdlibm.h"
#include <stdio.h>
int main(){
printf("%.16g\n" , sin( M_PI_4 ));
}
vbox$ cc -g zmain.c s_sin.o k_sin.o k_rem_pio2.o k_cos.o s_floor.o s_scalbn.o e_rem_pio2.o
関係者を招聘してコンパイル & gdb
Breakpoint 1, sin (x=0.78539816339744828) at s_sin.c:54 54 double y[2],z=0.0; (gdb) n 58 ix = __HI(x); (gdb) 61 ix &= 0x7fffffff; (gdb) 62 if(ix <= 0x3fe921fb) return __kernel_sin(x,z,0); (gdb) s __kernel_sin (x=0.78539816339744828, y=0, iy=0) at k_sin.c:66 66 ix = __HI(x)&0x7fffffff; /* high word of x */ (gdb) 67 if(ix<0x3e400000) /* |x| < 2**-27 */ (gdb) 69 z = x*x; (gdb) n 70 v = z*x; (gdb) 71 r = S2+z*(S3+z*(S4+z*(S5+z*S6))); (gdb) 72 if(iy==0) return x+v*(S1+z*r); (gdb) 74 } (gdb) sin (x=0.78539816339744828) at s_sin.c:78 78 } (gdb) s _libc_printf (fmt=0x1514c45d "%.16g\n") at /usr/src/lib/libc/stdio/printf.c:43 43 va_start(ap, fmt); (gdb) c Continuing. 0.7071067811865475
これが正確な値かな。
vbox$ bc scale=30 1 / sqrt(2) .707106781186547524400844362105
Algorithm
k_sin.c にアルゴリズムの説明が載ってた。
/* __kernel_sin( x, y, iy) * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854 * Input x is assumed to be bounded by ~pi/4 in magnitude. * Input y is the tail of x. * Input iy indicates whether y is 0. (if iy=0, y assume to be 0). * * Algorithm * 1. Since sin(-x) = -sin(x), we need only to consider positive x. * 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0. * 3. sin(x) is approximated by a polynomial of degree 13 on * [0,pi/4] * 3 13 * sin(x) ~ x + S1*x + ... + S6*x double S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
この係数って、テーラー展開のあれかな？
(%i2) taylor(sin(x), x, 0, 13);
3 5 7 9 11 13
x x x x x x
(%o2)/T/ x - -- + --- - ---- + ------ - -------- + ---------- + . . .
6 120 5040 362880 39916800 6227020800
x3の時の係数はS1となってる。テーラー展開のそれらを当てはめると、
1 / 6 = 0.166666666666666666666666666667 1 / 120 = 8.33333333333333333333333333333e-3 1 / 5040 = 1.98412698412698412698412698413e-4 1 / 362880 = 2.7557319223985890652557319224e-6 1 / 39916800 = 2.50521083854417187750521083854e-8 1 / 6226020800 = 1.60590438368216145993923771701e-10
こんな計算値になった(こういうの emacs/calcは便利)。微妙に数値が異なる のは、計算方法の違いを吸収する為の調整かな > 教えて Sunの人。って、も う、誰もいない。
pack/unpack by ruby
各定数に、コメントが記載されてた。一応確認する。
static double
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
int main(){ }
do gdb
Reading symbols from a.out...
(gdb) b main
Breakpoint 1 at 0x11b2: file dd.c, line 10.
(gdb) r
Starting program: /tmp/z/a.out
Breakpoint 1, main () at dd.c:10
10 int main(){}
(gdb) p S6
$1 = 1.5896909952115501e-10
(gdb) p/x S6
$2 = 0x3de5d93a5acfd57c
(gdb) x/2xw &S6
0x404050 <S6>: 0x5acfd57c 0x3de5d93a
リトルエンディアンな石です。
LLでやるには、どうする？
[sakae@arch z]$ ri
>> Array.pack
:
Float | Array |
Directive | Element | Meaning
---------------------------------------------------------------------------
D d | Float | double-precision, native format
F f | Float | single-precision, native format
E | Float | double-precision, little-endian byte order
e | Float | single-precision, little-endian byte order
G | Float | double-precision, network (big-endian) byte order
g | Float | single-precision, network (big-endian) byte order
まずはパックするんだな。そうするとバイナリーになる。そのバイナリーを16 進数に変換すればOK。
irb(main):001:0> v = 1.58969099521155010221e-10
=> 1.58969099521155e-10
irb(main):006:0> [v].pack("E").unpack("H*")
=> ["7cd5cf5a3ad9e53d"]
これって、インテル石の、いやがらせだな。
irb(main):007:0> [v].pack("E").reverse.unpack("H*")
=> ["3de5d93a5acfd57c"]
ならば反転してしまえ。
String.unpack では、
String | | Directive | Returns | Meaning ----------------------------------------------------------------- H | String | hex string (high nibble first) h | String | hex string (low nibble first) :
こんな仕掛けも用意されてた。日本的な、おもてなし。
対抗馬は、python
etc
to_lisp(); <=> (to-maxima) で、行ったりきたり。
rlwrap maxima すると、矢印キーで、履歴をたどって入力できるんで便利。