22 exflib を使ってみる (2018年秋)

(しばらく工事中)

以前のメモは

「多倍長計算ノート」

の第2節にある。

久しぶりに exflib を使ってみる。 動機は、これまた久しぶりに MATLAB を使うために色々見直していて、 それならば前から気になっていた MATLAB から exflib を使う、 というのを試してみよう、ということである。

まず g++, gfortran 向けの exflib

(MATLAB 用とは独立なので、MATLAB で使うだけならば、必要はないけれど。)

Mac用の最新版は、

http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib/exflib-x64-bin-20180620.tar.bz2

あ、以前と違って .pkg ではないのか。

curl -O http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib/exflib-x64-bin-20180620.tar.bz2
tar xjf exflib-x64-bin-20180620.tar.bz2
cd exflib-x64-bin-20180620
ls
less README.txt

sample-cxx, sample-f90, sample-par というディレクトリィが現れる。これ自体は以前と同じだ。

cd sample-cxx
make

あれ、ひっかっかった。 function というプログラムのコンパイルで、 exflib::work なんて知らないと言われる。 サンプル・プログラムの作成日時が古い。 作者本人は現在のバージョンの exflib で試していないのだと思われる。 それは抜けば make 出来る。

(Makefile をちょっと修正)

% diff Makefile Makefile.org
33,34c33
< #PROG=fact napier pi zero function test complex taylor contfrac matrix version subnormal
< PROG=fact napier pi zero test complex taylor contfrac matrix version subnormal
---
> PROG=fact napier pi zero function test complex taylor contfrac matrix version subnormal
% make

(無事終了)

出来上がった実行形式を起動すると、 一応動いているみたいだけれど、 どうなったら正しく動いているのか、すぐ分かるようになっていない。 「ソースを読め」ということらしい。 作者が動かした環境での結果くらいを含めておいてくれると楽なのだけど。

ライブラリィ・アーカイブ・ファイルとインクルード・ファイルを システムのディレクトリィにコピーして、 それを使ってコンパイルしてみる。

sudo cp -p ../exflib-x86_64-intelmac/libexfloat.a /usr/local/lib
sudo cp -p ../exflib-x86_64-intelmac/*.h /usr/local/include
g++ -Wall -O2 -m64 -I/usr/local/include -o napier napier.cpp -lexfloat
./napier

これで

2.718282e+0
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427

となった。まあ、いいのかな。-m64 が本当に必要かどうか分からない。 要らないような気がする。 少なくとも napier.cpp の場合、 つけてもつけなくても出来上がる実行形式に差はない。

exflib を使うプログラムをコンパイルするには、 次の2点に注意する。

(a)

exfloat.h, exfutil.h が インクルード出来るようにする (-I /usr/local/include などが必要になるかも)。

(b)

libexfloat.a をリンク出来るようにする。 (-lexfloat は必須。 もしかすると -L/usr/local/lib のように、 ディレクトリィを指定する必要があるかも。)

まとめると

こうするとインストールできるはず

curl -O http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib/exflib-x64-bin-20180620.tar.bz2
tar xjf exflib-x64-bin-20180620.tar.bz2
sudo cp -p exflib-x64-bin-20180620/exflib-x86_64-intelmac/libexfloat.a /usr/local/lib
sudo cp -p exflib-x64-bin-20180620/exflib-x86_64-intelmac/*.h /usr/local/include

プログラムの書き方については、 ずっと以前 「多倍長計算ノート」 の第2節に少し書いておいた。

では MATLAB 向けの exflib

http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib/exflib-MATLAB-20180902-distribute.zip

ターミナルで

curl -O http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib/exflib-MATLAB-20180902-distribute.zip
unzip exflib-MATLAB-20180902-distribute.zip
cd exflib-MATLAB-20180902-distribute
ls
cat README.txt
cd ..

MATLAB を起動して、実行してみる。

MATLABで

>> cd ${SOMEWHERE}/exflib-MATLAB-20180902-distribute
>> ls
>> edit fact
>> fact

これも一体どうなるべきものなのか、説明がない。

使うたびに cd したり、PATH を気にしたくないので、 設定しようと考えたのだけど、どうすべきか、説明が書かれてない。 自分でやんなさい、ということかな。

動くための条件は次の2つらしい。

(a)

exflib-MATLAB-20180902-distribute の下にある exflib の下にあるものを検索出来ること (addpath() を使って、 exflib を path に加えるのが自然)

(b)

exflib-MATLAB-20180902-distribute の下にあるメソッド・ディレ クトリィ @exfloat が検索できること。 メソッド・ディレクトリィは addpath() 出来ないので、 検索されるディレクトリィに移す必要がある。

結局次のようにした。 まず、

cd ~/Documents/MATLAB
cp -pr ${SOMEWHERE}/exflib-MATLAB-20180902-distribute/exflib .
cp -pr ${SOMEWHERE}/exflib-MATLAB-20180902-distribute/@exfloat exflib

のようにコピーしてから (MATLAB の下にexflib をコピーし、 その MATLAB/exflib の下に@exfloat をコピーした)、 startup.m に

addpath('~/Documents/MATLAB/exflib');

という設定を書き加えた。

以上をまとめると

こうするとインストール＆設定出来るはず (その1)

curl -O http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib/exflib-MATLAB-20180902-distribute.zip
unzip exflib-MATLAB-20180902-distribute.zip
cp -pr exflib-MATLAB-20180902-distribute/exflib ~/Documents/MATLAB
cp -pr exflib-MATLAB-20180902-distribute/@exfloat ~/Documents/MATLAB/exflib
echo "addpath('~/Documents/MATLAB/exflib')">> ~/Documents/MATLAB/startup.m

実は exflib-MATLAB-20180902-distribute の下に、 いくつか .m ファイルがあるので、 それも見えるようにしようと考えると、

こうするとインストール＆設定出来るはず (その2)

curl -O http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib/exflib-MATLAB-20180902-distribute.zip
unzip exflib-MATLAB-20180902-distribute.zip
mv exflib-MATLAB-20180902-distribute ~/Documents/MATLAB/exflib-20180902
echo "addpath('~/Documents/MATLAB/exflib-20180902')">>~/Documents/MATLAB/startup.m
echo "addpath('~/Documents/MATLAB/exflib-20180902/exflib')">>~/Documents/MATLAB/startup.m

ともあれ、テスト。

まず $100!$ を計算してみよう。

MATLABのコマンド・ウィンドウで

f=exfloat(1);
for n=1:100
  f=f*n;
  fprintf('fact %d = %s\n', n, num2str(f, '%f'));
end

精度は動的には変化させられず、@exfloat/exfloat.m で指定する。

exfloat.m の最初の8行

classdef exfloat

  properties (Constant, Access=private)

    %----------------------------------------
    % required precision in decimal digits
    %----------------------------------------
    precision10 = 500;

(応用数理学会でのデモの際は precision10 = 100; だったけれど、 変えてある。)

線型差分方程式 (符号に注意²)

$$a_n=\frac{34}{11}a_{n-1}-\frac{3}{11}a_{n-2}$$

の一般解は

$$a_n=\frac{A}{11^{n-1}}+B\cdot 3^{n-1}$$   ($A$, $B$ は任意定数)

であり、初期条件

$$a_1=1,\quad a_2=\frac{1}{11}$$

を満たす解は $A=1$, $B=0$ として得られる $a_n=\dfrac{1}{11^{n-1}}$ である。

clear classes;
a(1)=exfloat(1);
a(2)=exfloat('1/11');
for n=3:50
  a(n)=exfloat('34/11')*a(n-1)-exfloat('3/11')*a(n-2);
end
plot(a)

figure
semilogy(a)

exfloat の代わりに普通の double では

clear('a')
a(1)=1;
a(2)=1/11;
for n=3:50
  a(n)=34/11 * a(n-1)-3/11*a(n-2);
end
plot(a)
figure
semilogy(a)

ここでは片側対数目盛りのグラフのみ載せておく。

図: exfloat では

図: double では

$a_n=\dfrac{1}{11^{n-1}}$ であるから、 まっすぐ右下がりの直線となるべきものである。 double で計算した場合は、 $n=13$あたりから丸め誤差の影響が現れている。