more maxima
東大の入試問題。
円周率が、3.05 より大きい事を証明せよ
まるで、電車の中吊り広告みたい。オイラーは降参です。解説を見ると、
半径１の円に内接する正６角形は、その辺の長さがおのおの１だから、全周は６。直径の２で割れば、パイは３。では、正８角形はどうよ？
頂角が45度の扇形になる。ピザやケーキを８等分して、円弧の部分を切り落とし、その底辺の長さが分かれば、それを８倍して全周が求まる。
頂点から、垂線を下ろすと、角度が22.5度の直角三角形が２つ出来る。その斜辺は半径の１に等しい。よって、sin(22.5度)を求めて、それを２倍すれば、８角形の一辺の長さが求まる。
試験会場にはsinの表も電卓も無いんで、何とか頭をひねれ。三角関係の半角の法則があるそうな。
2 a 1 - cos a
sin --- = ------------
2 2
cos(45度)は、1/sqrt(2) = sqrt(2)/2 になる事ぐらいは受験生の常識。これで、何とか式の変形と筆算で行けそう。オイラーは受験生じゃないんで、楽します。
(%i4) 8 * 2 * sin(%pi / 8);
%pi
(%o4) 16*sin(---)
8
(%i5) float( % / 2);
(%o5) 3.061467458920718
maximaが有れば、東大もスイスイ？
受験生は、16 * sin(22.5) = 16 * sqrt( (1 - cos(45) / 2) のように、左式を開平。同様に右式も開平。cos(45)に、sqrt(2)/2を適用。こうしておいて計算する訳だ。大変なこっちゃなあ。
難問
a^3 + b^3 + c^3 = 42
上に式を満たす整数 a,b,cを求めよ。（答えは巻末に）
Bug発見
前回の円周率を求めるやつに、うっかりBugが有った。
mypi(A,B) := (soto(A,B) + naka(A,B)) / 4$ to mypi(A,B) := (A + B) / 4$
一回余分にやってたので、訂正します。訂正後のso.mac実行の冒頭部分。以前は、3.152...から結果が表示されてた。
3.187587978952741
3.152021515166292
3.144136698987599
:
同様にmy.macの方も修正。
(%i4) float( my(8, 4*sqrt(2),1) ); (%o4) 3.187587978952741
これが、正８角形でのパイの予想。
(%i6) float( my(8, 4*sqrt(2),11) ); Maxima encountered a Lisp error: arithmetic error FLOATING-POINT-OVERFLOW signalled Automatically continuing. To enable the Lisp debugger set *debugger-hook* to nil.
余り欲張ると、エラーになる。
(%i7) float( my(8, 4*sqrt(2), 20) ); Heap exhausted during garbage collection: 0 bytes available, 16 requested. Gen Boxed Unboxed LgBox LgUnbox Pin Alloc Waste Trig WP GCs Mem-age fatal error encountered in SBCL pid 5889(tid 0x7fd677df5b80): GC invariant lost, file "gencgc.c", line 472 Welcome to LDB, a low-level debugger for the Lisp runtime environment.
更に欲張るとmaximaが検出するエラーを通り越して、地のsbclのマシン語レベルのdebuggerが頭をもたげて来た。地獄に落ちて閻魔様がお出ましになった状態だな。メモリー欠乏症って言う根源に触れちゃったんだ。
ldb> backtrace Backtrace: 0: COPY-TREE, pc = 0x5293866a, fp = 0x7fd676bb4338 1: COPY-TREE, pc = 0x52938631, fp = 0x7fd676bb4368 : 24: COPY-TREE, pc = 0x52938631, fp = 0x7fd676bb47b8 25: MAXIMA::SIMPTIMES, pc = 0x5212f28f, fp = 0x7fd676bb4820 26: MAXIMA::MEVALARGS, pc = 0x521a5246, fp = 0x7fd676bb4878 : 96: MAXIMA::MEVAL1, pc = 0x5216e49d, fp = 0x7fd676bb5e48 97: MAXIMA::MEVAL, pc = 0x5210dca4, fp = 0x7fd676bb5e68 98: MAXIMA::MLAMBDA, pc = 0x52224a73, fp = 0x7fd676bb5f08
そもそも、浮動小数の演算でオーバーフローさせちゃうようなコードが悪いんだろう。 myが呼び出される度に、数式が肥大化していくからねえ。肥大化しない対策を打とう。
soto(A,B) := 2 * A * B / (A + B)$
naka(A,B) := sqrt(soto(A,B) * B)$
mypi(A,B) := float( (A + B) / 4)$
my(A,B,n) := block(
if n = 0 then
return( mypi(A,B) )
else
return( my(float(soto(A,B)), float(naka(A,B)), n - 1) )
)$
再帰する時、引数を数値にしちゃうと言う対策を施しました。まあ、こうしちゃうと、誤差が累積する事は目に見えているんですけどね。
(%i7) my(8, 4*sqrt(2), 20); (%o7) 3.141592653589941 (%i8) my(8, 4*sqrt(2), 30); (%o8) 3.141592653589793 (%i9) float(%pi); (%o9) 3.141592653589793 (%i12) %o9 - %o7; (%o12) - 1.474376176702208e-13
本格的な競技的パイ競争でもしない限り、20回の試行でも十分な精度が得られるな。
多倍長浮動小数点（ビッグフロート）
円周率とかを多数桁求めようとした時、多分必要になるであろう機能。 通常の浮動小数表現は16桁程度。多倍長精度浮動小数点がサポートされてる。使う前に、fpprecと言う変数に、仮数の桁数を指定しておく。
(%i13) 1.0/7; (%o13) 0.1428571428571428 (%i15) fpprec:50; (%o15) 50 (%i16) 1.0b0/7; (%o16) 1.4285714285714285714285714285714285714285714285714b-1
ビッグフロートを指定する時は、仮数b指数 のように指定する。良く 1.23e2 とかやる奴の仮数部拡張版って事だ。
(%i20) bfloat(%pi);
(%o20) 3.1415926535897932384626433832795028841971693993751b0
(%i21) rationalize(%);
587704679302172021937255065567143824442131843125525
(%o21) ---------------------------------------------------
187072209578355573530071658587684226515959365500928
浮動小数を、有理数に変換なんて事も出来る。
いろいろな円周率をビッグフロートにするのもいいな。
便利な奴
maximaから出力される整形された式を、maximaが認識される入力形式に変換したい事がある。
(%i5) solve(2 / a = 1 / A + 1 / B, a), A=8, B=4*sqrt(2);
13/2
2
(%o5) [a = --------]
5/2
2 + 8
(%i6) grind(%);
[a = 2^(13/2)/(2^(5/2)+8)]$
(%o6) done
(%i7) [a = 2^(13/2)/(2^(5/2)+8)];
13/2
2
(%o7) [a = --------]
5/2
2 + 8
上で、%o5がそれだどしよう。grind関数に渡してあげるだけで、入力形式が得られる。 なんか、こんなのschemeにもあったな。
(%i8) ibase: 16; (%o8) 16 (%i9) ff; (%o9) ff (%i10) 0ff; (%o10) 255
任意な進数の設定。16進の場合は頭に0を補って、数字のふりをする必要がある（場合があります） ibaseが有るなら、obaseもあるだろう。勿論ありますよ。
many Common Lisp
例によってmaximaのdirを漁っていたんだ。そしたら面白い物を発見。
from maxima/README.lisps
Comparison of execution times (in seconds) for the run_testsuite() function for Maxima 5.36.0 as reported on the Maxima mailing list: 64 Bit (Gentoo Linux): gcl-2.6.12 152 sbcl-1.2.10 155 ccl-1.10 313 ecl-15.3.7 559 clisp-2.49 1060 32 Bit (Gentoo Linux) sbcl-1.2.10 170 gcl-2.6.12 177 cmucl-20e 253 ccl-1.10 293 ecl-15.3.7 556 clisp-2.49 844
以前sbclを作る時clispが必要って事で嘆いたけど、ここでも嘆きですよ。それはさておき、gclがsbclと遜色無いスピードで動く事に感心。ベンチマークの為のコードだと、チューニングして、世間を欺く事も出来ようが、今回はそれもなさそう。
CommonLispと言えども、微妙な違いがあるだろうから、sbcl一択から脱して、gclも試してみるかな。 A Road to Common Lisp (Common Lisp への道) こういう指南書も有る事だしね。
gcl
debian機で試す。 GCL - GNU Common Lisp ソースを頂いてきた。既に開発は終了しちゃってるのかな？
で、maxima用にgclを作るには、config時に、 --enable-ansi しておけとの事。KCLが昔風のやつで、それを世界標準にグレードアップするフラグなんだな。ああ、ANSIってアメリカ版のJIS規格か。
で、コンパイルエラーですよ。正確にはリンクエラー。エラー潰しの元気は無いので、あらぬ方向に走る事にします。
を参考にソースからbuild。cd src して、その中でapt source gcl
-rw-r--r-- 1 root root 7031837 Oct 29 2014 gcl_2.6.12.orig.tar.gz -rw-r--r-- 1 root root 1832 Feb 6 2019 gcl_2.6.12-83.dsc -rw-r--r-- 1 root root 416484 Feb 6 2019 gcl_2.6.12-83.debian.tar.xz
こんなやつをsrc下にぶちまけてから、debian用にパッチの当たった、 gcl-2.6.12/が出来た。その中に移動してbuild。
また、外側に、結果のファイルが、 gcl_2.6.12-83_i386.deb gcl_2.6.12-83_i386.{build,buildinfo,changes} な形で出来上がった。後は、インストール。
起動したら、古い版で動いた。はてなと思って、/usr/bin/gclを覗くと
if [ "$GCL_ANSI" = "" ] ; then
EXE=saved_gcl;
else
EXE=saved_ansi_gcl;
fi
こんな事が書いてある。根本で設定するには、/etc/default/gclに書いておけとな。
The Debian package gcl ---------------------- GCL is one of the oldest free common lisp systems still in use. Several production systems have used it for over a decade. The common lisp standard in effect when GCL was first released is known as "Common Lisp, the Language" (CLtL1) after a book by Steele of the same name providing this specification. Subsequently, a much expanded standard was adopted by the American National Standards Institute (ANSI), which is still considered the definitive common lisp language specification to this day. Debian GCL now installs both the small 'traditional' lisp image designed to conform to a pre-ANSI Lisp standard, and an experimental ANSI image. Please note that ANSI support in GCL is still preliminary. On an ansi-test suite written by a GCL developer, GCL fails on a little under 3 percent of the tests. Details can be found in /usr/share/doc/gcl/test_results.gz. To toggle the use of the ANSI image, set the environment variable GCL_ANSI to any non-empty string.
debianのgcl releaseにも書いてあった。
これでmaximaをコンパイル出来るな。面倒なのでconfigure方式で作成した。で、試運転。 dumpされた物は、gcl + maximaの塊になってるので、直接起動。maximaを起動するスクリプトを見ると、GCがドーダコーダだから、環境変数で細かく制御しなければいけないようだけど、超簡略版での起動。
debian:~$ cd src/maxima-5.43.0/src/
debian:src$ binary-gcl/maxima
GCL (GNU Common Lisp) 2.6.12 ANSI Fri Apr 22 15:51:11 UTC 2016
Source License: LGPL(gcl,gmp), GPL(unexec,bfd,xgcl)
:
>(run)
Maxima 5.43.0 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.12
:
(%i1) sin(%pi/4);
1
(%o1) -------
sqrt(2)
例によって(run)でmaximaが動き出す。
(%i2) :lisp (trace simp-%sin %piargs-sin/cos)
(SIMP-%SIN %PIARGS-SIN/COS)
(%i2) sin(%pi/4);
1> (SIMP-%SIN ((%SIN) ((MTIMES SIMP) ((RAT SIMP) 1 4) $%PI)) 1 NIL)
2> (%PIARGS-SIN/COS ((MTIMES SIMP) ((RAT SIMP) 1 4) $%PI))
<2 (%PIARGS-SIN/COS ((MEXPT SIMP) 2 ((RAT SIMP) -1 2)))
<1 (SIMP-%SIN ((MEXPT SIMP) 2 ((RAT SIMP) -1 2)))
1
(%o2) -------
sqrt(2)
sbclの時と違ってスッキリ感がする。その代わりシンボル名が大文字になってるのは、ちょいと馴染めないな。昔はこうだったんだね。
そんじゃ、slimeをgclで動かして、、、
>;; Loading "/home/sakae/.emacs.d/elpa/slime-20190930.1713/swank-loader.lisp"
Error:
Fast links are on: do (si::use-fast-links nil) for debugging
Signalled by FUNCALL.
SIMPLE-ERROR: Too few arguments are supplied to defmacro-lambda-list.
Broken at FUNCALL. Type :H for Help.
1 Return to top level.
SWANK-LOADER>>
見事にエラー。調べたら、 GCL and slimeだれも、gclなんて動かさないよー。機能不足で誰もしり込みしてるとの事。 slimeの本山でも、サポート対象にはなっていない。残念ではありますがね。
何とかslime無しで動かんか？
真似をしてみる。
(%i2) to_lisp(); Type (to-maxima) to restart, ($quit) to quit Maxima. MAXIMA> (compile-file "trigi.lisp") #p"/home/sakae/src/maxima-5.43.0/src/trigi.o" MAXIMA> (load #p"/home/sakae/src/maxima-5.43.0/src/trigi.o") 82019
ソースファイルを別の端末で修正して、それをコンパイルしてからロード。特定の関数だけコンパイルしたら、リンクエラーだった為。でも、走らせると、maximaのトップレベルで落ちちゃった。相性悪いんすかね？
gclをコンパイルしたdirの中にelispってのが有り、maximaと連携出来る風には、説明されてるんだけどね。折角gclで動く物を作ったけど、slimeで動かんのは致命的。
悔し紛れに、ベンチマークを取って、終わりにするか。
(%i1) run_testsuite();
Testsuite run for GNU Common Lisp (GCL) GCL 2.6.12:
Running tests in rtest_rules: 119/119 tests passed
Running tests in rtestnset: 617/617 tests passed
:
No unexpected errors found out of 11,637 tests.
Tests that were expected to fail but passed:
/tmp/maxima-5.43.0/tests/rtest_gamma.mac problem:
(307)
real time : 290.959 secs
run-gbc time : 262.769 secs
child run time : 0.930 secs
gbc time : 22.370 secs
(%o0) done
今度はsbcl環境。
SBCL 1.4.16.debian + Maxima 5.43.0 : No unexpected errors found out of 11,637 tests. Evaluation took: 215.232 seconds of real time 213.235529 seconds of total run time (211.465734 user, 1.769795 system) [ Run times consist of 5.504 seconds GC time, and 207.732 seconds non-GC time. ] 99.07% CPU 23,330 forms interpreted 15,154 lambdas converted 472,344,619,791 processor cycles 22,369,824,584 bytes consed (%o0) done
maximaのREADMEでは、ほぼ互角って事だったけど、290秒 対 215秒の差が付いたな。
参考までに、CentOS上で SBCL 1.5.6 + Maxima 5.43.0 要する64Bit環境ね。CentOSには、2.5Gを割り当てている。メモリーを増やすとGCが減るのかな？
No unexpected errors found out of 11,637 tests. Evaluation took: 85.551 seconds of real time 84.309844 seconds of total run time (82.722969 user, 1.586875 system) [ Run times consist of 3.227 seconds GC time, and 81.083 seconds non-GC time. ] 98.55% CPU 21,196 forms interpreted 16,194 lambdas converted 205,320,109,930 processor cycles 4 page faults 35,340,793,888 bytes consed
slime
slimeではgclをサポートしてなくて残念。もう一つ残念なのはschemeをサポートしていない事。昔、一生懸命に探した事があったけど、あれから有志は現れたのだろうか？
What is the closest thing to Slime for Scheme?
どうも無理っぽい。新版のlispであるclojureは、ciderってのに移行しちゃったし、やっぱりschemeはメガヒットが無いから、slimeする動機にならないんだろうね。
Common Lisp の開発環境 slime こちらさんは、べったりCLだしね。
apt install cl-swank とかして /usr/share/common-lisp/source/slime/start-swank.lisp を動かせばよい。clisp でも sbcl でもどちらでもよい。 というか通常は自動起動するから(vim や emacs から) そもそも起動しなくてもよい。
viな人もslimeを使おうとされてるし。
elp/slime-2019xxx
contrib/slime-scheme.el ;;; slime-scheme.el --- Support Scheme programs running under Common Lisp ob$ ls contrib/*scm swank-kawa.scm swank-mit-scheme.scm swank-larceny.scm swank-r6rs.scm
こういう物も有るんだ。知らんかった。
回答
(%i4) a:-80538738812075974; (%o4) - 80538738812075974 (%i5) b:80435758145817515; (%o5) 80435758145817515 (%i6) c:12602123297335631; (%o6) 12602123297335631 (%i7) a^3+b^3+c^3; (%o7) 42
検算した人居なかったので、maxima君にさせたよ。
それにしても、17桁の数字、万億兆京の位とは、恐れ入る。日本の国家予算が100兆円台、そして借金が1000兆円って事からすれば、いかに巨大な数か。
maxima君は、こんな有り得ないような数字の計算をいとも容易くやってのける。素敵な奴だなあ。