maximaの使用例
絵本作家の かこさとしさんの『科学者の目』なんて本を読んだんだ。 41人の科学者が紹介されてた。どんな視点で科学に臨んだかを簡潔に紹介している。 アインシュタインとかダ・ビンチとか、日本からだと長岡半太郎先生とか。
女性では、丹下ウメさん。新聞に連載中に中学生から、日本にキュリー夫人みたいな素敵な科学者はおられますかって投書があった。それに答えて、かこさんが推薦されたそうな。
ノーベル賞のノーベルさんも出てた。物理、化学、医学、経済、平和賞。それに異端児のように文学賞が加わっている。青年の頃の彼は文学に入れ込んでいたので、その名残だそうだ。知らんかった。
それにしても、科学系中心の賞なのに、科学の母である数学が無いのは何故？ そこまでの説明は、子供向けと言う事で省いたのかな。 聞く所によると、数学を憎んでいたからーー。これ有名な話みたいですね。チコちゃんが聞いてくるぞ。必須雑学ですよ。ノーベル数学賞はなぜ存在しない？
中国の昔(429-500)の科学者、祖沖之 が、パイを正確に計算してたって事も紹介されてた。
アルキメデスもパイ好きで、22/7と223/71の間に挟まっているよと計算してる。祖さんは、それとは独立に情熱を傾けて求めている。
22/7を約率、355/113を密率として、導きだした。
導出方法は、円に外接する多角形（の辺の合計）を、Aとする。内接する多角形の方をBとする。 同様にその多角形の2倍の多角形のそれを、それぞれ、a,bとすると、
A > 2πr > B 2/a = 1/A + 1/B b = sqrt(aB)
という関係が成り立つ。これを使って、倍々に多角形の角数を増やしていった。密数に達するまでの精度を得るには、192角形が必要とか。
結果をmaximaで確認。
(%i6) float(22/7); (%o6) 3.142857142857143 (%i7) float(355/113); (%o7) 3.141592920353983 (%i8) float(%pi); (%o8) 3.141592653589793
355/113は、科学計算に使うには十分過ぎる精度を持っている。パイへの情熱恐るべし。
debugger for maxima
maximaのマニュアルを見ていたんだ。そしたら上級者が使って便利なdebuggerが紹介されてた。maxima語でスクリプトを書くような人の最終兵器だからね。
早速、概要を眺める。って言っても、プロンプトから:hとかするだけだけだけど。
(%i1) :h
Break commands start with ':'. Any unique substring may be used,
eg :r :re :res all work for :resume.
Command Description
----------- --------------------------------------
:break Set a breakpoint in the specified FUNCTION at the
specified LINE offset from the beginning of the function.
:
:lisp Evaluate the lisp form following on the line
:quit Quit this level
この陣容を見て、ちょいとしたアイデアが沸いてきた。
(%i1) :lisp (trace simp-%sin)
(SIMP-%SIN)
(%i1) sin(%pi/4);
0: (MAXIMA::SIMP-%SIN ((MAXIMA::%SIN) ((MAXIMA::MTIMES MAXIMA::SIMP) ((MAXIMA::RAT MAXIMA::SIMP) 1 4) MAXIMA::$%PI)) 1 NIL)
0: SIMP-%SIN returned ((MEXPT SIMP) 2 ((RAT SIMP) -1 2))
1
(%o1) -------
sqrt(2)
(%i2) :lisp (untrace simp-%sin)
T
(%i2) sin(%pi/4);
1
(%o2) -------
sqrt(2)
trace/untraceの設定が、maximaのプロンプト上から、超簡単に出来る。これは泣いて喜ぶべき事ですよ。rmaxima上から使えば、コマンド履歴を辿って、修正実行も極めて簡単に出来るからね。
slime-edit-definition (M-.)
slimeを起動しておいて、C-c C-l zzz.lispして、maximaのコードをローディングする。それから、下記のように、packageをmaximaに切り替える。
CL-USER> (in-package maxima) #<PACKAGE "MAXIMA"> MAXIMA>
下記は、trigo.lispのtailの部分。（別にどんなファイルの何処でもいいんだけど）
(defun mpc1 (dl ul result f di ui) (cond ((= 0 ui) (cons (reverse dl) result)) :
やおら、reverseにカーソルを合わせてから、M-. して、コードの定義元を探してみる。
/home/christophe/src/sbcl-release-dir-20190826-lVOgdbu6p/sbcl-1.5.6/src/code/list.lisp (DEFUN REVAPPEND) /home/christophe/src/sbcl-release-dir-20190826-lVOgdbu6p/sbcl-1.5.6/src/compiler/fndb.lisp (DECLAIM REVAPPEND SB-C:DEFKNOWN)
画面が割れて、上記のような案内が出てきた。かまわずn,pで先を辿ると
Use M-x make-directory RET RET to create the directory and its parents
slimeが悲鳴を上げてきた。苦し紛れにdirを作るかですって。想像出来るように、christopheさんの作業場所だ。そこのsrc以下でsbcl本体をコンパイルしたんだな。残念ながら、ソースは見られない。何とかしたいぞ。
その前に、lispのコンパイル結果に、出目が刻印されてるのかな？ ちょっと探ってみる。
(base) [sakae@c8 binary-sbcl]$ pwd /usr/local/lib/maxima/5.43.0/binary-sbcl (base) [sakae@c8 binary-sbcl]$ strings maxima.core | grep christophe shin-christophe-2019-08-26-10-02-18 /home/christophe/src/sbcl-release-dir-20190826-lVOgdbu6p/sbcl-1.5.6/obj/asdf-cache/sb-bsd-sockets/constants.lisp-temp :
やっぱり、埋め込んであった。coreになる前の trigo.faslでも
(base) [sakae@c8 binary-sbcl]$ pwd /home/sakae/src/maxima/src/binary-sbcl (base) [sakae@c8 binary-sbcl]$ strings trigo.fasl | grep sakae compiled from "/home/sakae/src/maxima/src/trigo.lisp" d%/home/sakae/src/maxima/src/trigo.lisp
埋め込んである。CLのアプリを特殊な用途に使う場合、出目が分からないようにしましょう。そうしないと、思わぬ所から足がつきまっせ！
install sbcl on OpenBSD 6.6
そんなこんなでsbclをコンパイルする事にした。どうせやるならOpenBSD 6.6 にしてからだな。何も難しい事は無い。ports/lang/sbclの下にあるメニューを消化するだけ。
Makefileの事前調査をしたら、コンパイルするにはsbclが必要。そんな鶏/卵問題を解決する為には、種になるsbclの代わりにclispが必要。clispかあと思っちゃったよ。一時的に入れておいて、用が済んだら消せばいいか。他に必要なのは、gmakeとtexinfoですって。これが入っていないと、勝手にこれを先にコンパイルするからなあ。事前にpkg_addで入れておいたよ。 そしたらtexinfoの付属品としてclispも付いてきた。
81610 p1 I+p 0:08.81 make 55826 p1 I+p 0:00.02 /bin/sh -ec lock=sbcl-1.5.5; export _LOCKS_HELD=" sbcl-1.5.5"; /usr/bin/perl /usr/ports/infrastructur 52083 p1 I+p 0:00.05 make _internal-all 21929 p1 I+p 0:00.02 make do-build 51160 p1 I+p 0:00.02 /bin/sh -ec cd /usr/ports/pobj/sbcl-1.5.5/sbcl-1.5.5 && /usr/bin/env -i GNUMAKE=gmake INSTALL_ROOT=/usr/ 44056 p1 I+p 0:00.01 /bin/sh make.sh --prefix=/usr/local --with-sb-core-compression --with-sb-xref-for-internals --xc-host=/u 8584 p1 R/3 29:11.42 /usr/local/lib/clisp/base/lisp.run -B /usr/local/lib/clisp -M /usr/local/lib/clisp/base/lispinit.mem -N 88324 p1 I+p 0:00.02 sh make-host-2.sh
のんびりとclispが走って、一生懸命にsbclを作っていますよ。
ob$ sbcl This is SBCL 1.5.5.openbsd, an implementation of ANSI Common Lisp. More information about SBCL is available at <http://www.sbcl.org/>.
出来上がったので、取り合えずの試走です。実技試験は、maximaのコンパイル。無事完了。
そして、CentOS 8に入ってるsbclでは出来なかった事を確認。コンパイルした跡地を放置したので、ちゃんとソースも温存されています。
/usr/ports/pobj/sbcl-1.5.5/sbcl-1.5.5/src/code/seq.lisp (DEFUN REVERSE) /usr/ports/pobj/sbcl-1.5.5/sbcl-1.5.5/src/compiler/seqtran.lisp (:DEFTRANSFORM REVERSE (LIST) "optimize") (:DEFTRANSFORM REVERSE (VECTOR) "optimize") /usr/ports/pobj/sbcl-1.5.5/sbcl-1.5.5/src/compiler/knownfun.lisp (:DEFOPTIMIZER REVERSE SB-C:DERIVE-TYPE) /usr/ports/pobj/sbcl-1.5.5/sbcl-1.5.5/src/compiler/fndb.lisp (DECLAIM REVERSE SB-C:DEFKNOWN)
それでは、ご対面ーーん。
;;;; REVERSE
(defun reverse (sequence)
"Return a new sequence containing the same elements but in reverse order."
(declare (explicit-check))
(seq-dispatch-checking sequence
(list-reverse sequence)
(vector-reverse sequence)
:
Smalltalkは、システムを作り上げているコードも容易に確認出来るそうだ。これを持って、同じ事が出来るようになったな。
それから、OpenBSD 6.6 にしたついでに、.emacs.d/init.el の澱を整理した。何とかモードは、極力削減。
;; slime (setq inferior-lisp-program "sbcl") (require 'slime) (slime-setup '(slime-repl slime-fancy slime-banner slime-company)) (add-hook 'slime-repl-mode-hook (lambda () (paredit-mode +1))) (add-hook 'slime-mode-hook (lambda () (paredit-mode +1)))
read-eval-print
ちょいと寄り道しちゃったけど、maximaの実行の流れを見てみる。今回は、continueから呼び出される、print系。前回作った美しい地図をじっと眺めて、予想をたてる。
displaが出力系の元締めで、その下に４種の系統がある。名前からして、出力が配列系とかlinear系とか、一般系とかに分かれている。配列系は、更にその下が５個に分類されてる。
今回はsinの流れを見るので、多分、普通系のoutputだろうと予測。traceで網を張る。
(%i2) :lisp (trace displa output)
(DISPLA OUTPUT)
(%i2) sin(%pi/4);
0: (MAXIMA::DISPLA ((MAXIMA::MLABEL) MAXIMA::$%O2 ((MAXIMA::MEXPT MAXIMA::SIMP) 2 ((MAXIMA::RAT MAXIMA::SIMP) -1 2))))
1: (MAXIMA::OUTPUT ((MAXIMA::D-HBAR 7) (-7 0) (-4 -1 #\) #\2 #\( #\t #\r #\q #\s) (3 1 #\1) (30 0) #\ #\) #\2 #\o #\% #\() 0)
1
(%o2) -------
sqrt(2)
1: OUTPUT returned NIL
0: DISPLA returned NIL
予想通りに、通過しましたね。これで、見る所が限定された。
displaを探すとdispla.lispに見つかった。
(defun displa (form &aux #+kcl(form form))
(if (not #.ttyoff)
:
(unwind-protect
(progn
(setq form (dimension form nil 'mparen 'mparen 0 0))
(checkbreak form width)
(output form (if (and (not $leftjust) (= 2 lines))
(- linel (- width bkptout))
0)))
;; make sure the linearray gets cleared out.
(fill linearray nil))))))
unwind-protectで囲まているって事は、予期せぬ失敗に備えているんだな。更に、このdisplaにカーソルを合わせて、C-c C-w C-c し、呼び出し元を列挙すると、
MAXIMA-DISPLAY /home/sakae/src/maxima/src/macsys.lisp CONTINUE CONTINUE $BUG_REPORT-IMPL MBREAK-LOOP /home/sakae/src/maxima/src/comm.lisp $DISPTERMS-IMPL :
色々な所から呼び出されている。汎用ルーチンなんだな。ざっと見して念力で、この経路だろうと当たりを付ける。地道に、displaに(break)を置いて、経路を炙り出す。どちらでもお好きな方を選択してください。
example
そんじゃ、例として冒頭にあげた祖さんの軌跡を追ってみるか。（おっ）パイに食らいつけって訳ね。
きっと開成中学の入試にも出るに違いない。半径１の円がぴったり収まる正方形と、円の中にぴったり入る正方形を考え、円の長さを推測し、パイの値をもとめなさい。この方法は、江戸時代に関さんが考えました。
(%i1) A: 8;
(%o1) 8
(%i2) B: 4 * sqrt(2);
5/2
(%o2) 2
(%i4) (A + B) / 4;
5/2
2 + 8
(%o4) --------
4
(%i5) float(%o4);
(%o5) 3.414213562373096
外接する正方形の一辺の長さは２になる。４辺で８．内接する正方形は、２辺の長さが１の直角三角形が４つ隣合わせになってると見做せる。直角三角形の斜辺の長さは、アルキメデスの法則により、sqrt(2)となる。
円周は、外接正方形と内接正方形と半分と予想される。それを２で割ればパイになるはず。
それじゃ、AとBを元に、正８角形の外周を求めてみる。
(%i8) solve(2 / a = 1 / A + 1 / B, a);
13/2
2
(%o8) [a = --------]
5/2
2 + 8
(%i9) float(%o8);
(%o9) [a = 6.627416997969525]
(%i10) b: sqrt(%o9 * B);
5/4 5/4
(%o10) [2 sqrt(a) = 2.574377011622331 2 ]
(%i15) float( (%o9 + b) / 4);
(%o15) [0.25 (a + 2.378414230005442 sqrt(a)) = 3.187587978952741]
3.41が3.18と段々近づいているな。以下、同じ事を繰り返せば精度が上がるはず。
%i8は、指示された方程式を、aについて解くための方法。事前に代入されてるAやBの値が使われて、簡略されてる。%i9は、前の結果%o8を使って数値化したものだ。
変数が多数ある式をそのまま解いて、後で各変数に具体的な数値を割り当てる方法もある。
(%i16) solve([a*x+b*y=p,c*x+d*y=q],[x,y]),a=1,b=3,c=2,d=4,p=23,q=34; (%o16) [[x = 5, y = 6]]
式とか求解したい変数が複数ある場合、鍵カッコで囲って指定する。連立方程式は、これでいけるよ。
(%i18) solve(2/x = 1/Y + 1/Z, x);
2 Y Z
(%o18) [x = -----]
Z + Y
変数に代入されてる値を利用しないで、記号的に解きたい時は、新しい記号を使うだけだ。 remvalue(A); とか、remvalue(A); して、削除していもいいけどね。
python方面には、記号で演算出来るsympyが有るけど、いちいち使う記号を宣言しないといけないので、非常に面倒だ。
(%i2) soto(A,B) := 2 * A * B / (A + B);
2 A B
(%o2) soto(A, B) := -----
A + B
(%i5) a: soto(8,4*sqrt(2));
13/2
2
(%o5) --------
5/2
2 + 8
(%i6) float(a);
(%o6) 6.627416997969525
(%i7) naka(A,B) := sqrt(soto(A,B) * B);
(%o7) naka(A, B) := sqrt(soto(A, B) B)
(%i8) b: naka(8, 4*sqrt(2));
9/2
2
(%o8) --------------
5/2
sqrt(2 + 8)
(%i9) (a + b) / 4;
13/2 9/2
2 2
-------- + --------------
5/2 5/2
2 + 8 sqrt(2 + 8)
(%o9) -------------------------
4
(%i10) float(%);
(%o10) 3.187587978952741
基本多角形の２種の辺長さA,Bを使って、倍多角形の辺長を求める関数、soto(A,B)とnaka(A,B)を定義してみた。ごちゃごちゃしてるのは、検算だ。
(base) [sakae@c8 tmp]$ cat so.mac soto(A,B) := 2 * A * B / (A + B)$ naka(A,B) := sqrt(soto(A,B) * B)$ mypi(A,B) := (soto(A,B) + naka(A,B)) / 4$ A : 8$ B : 4 * sqrt(2)$ for i:1 thru 10 do block( a : soto(A,B), b : naka(A,B), disp(float(mypi(a,b))), A : a, B : b)$
十分な精度まで達するかな。
(%i1) load("so.mac");
3.152021515166287
3.144136698987593
3.142224771100327
3.141750440437613
3.141632085156619
3.141602510535152
3.141595117767005
3.141593269630368
3.14159280759969
Maxima encountered a Lisp error:
arithmetic error FLOATING-POINT-OVERFLOW signalled
Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.
浮動小数もオーバーフロー（アンダーフロー）するんだ。 debian機ではどうかな。GCL環境だったはずだから。
debian:tmp$ maxima -b so.mac
Maxima 5.42.1 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.12
(%i1) batch("so.mac")
:
3.141593269630374
3.141592807599708
#<-nan>
最後は、計算出来ませんですって。理由無しで打ち切られるのは辛い。
で、上のforで回すやつを書いていて、微妙にイライラしてました。これって、格好の再帰関数の題材だよね。書いてみる。
(base) [sakae@c8 tmp]$ cat my.mac
soto(A,B) := 2 * A * B / (A + B)$
naka(A,B) := sqrt(soto(A,B) * B)$
mypi(A,B) := (soto(A,B) + naka(A,B)) / 4$
my(A,B,n) := block(
if n = 0 then
return( mypi(A,B) )
else
return( my(soto(A,B), naka(A,B), n - 1) )
)$
明示的にreturnを入れないと、だんまりの症状になるので注意。
(%i1) load("my.mac");
(%o1) my.mac
(%i2) float( my(8, 4 * sqrt(2), 1) );
(%o2) 3.152021515166287
(%i3) showtime: true;
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o3) true
(%i4) float( my(8, 4 * sqrt(2), 9) );
Evaluation took 1.1250 seconds (1.1370 elapsed) using 424.340 MB.
(%o4) 3.14159280759969
再帰関数からは、mypiの計算式が返って来るので、数値化してる。そのまま表示させちゃうと、大変（長い数式が表示されてしまい）やっかいな事になる。莫大な資源も使ってるしね。