maximaとか

一つ 頓智と言うかクイズをば。

有限なのに無限なものはなーーんだ？

それはー、図書館ですだ。蔵書は有限なのに、おいらが全部読もうと思ったって、生涯のうちに 全部は読めないから！ 今の所硬軟合せて、10冊/月 ぐらいのペースで読んでますけど。

先日町の図書館へ行ったら、ご自由にお持ち下さいの棚が用意されてた。図書を貰ってっても 良いって事？ いぶかしげに本を見ると、『除籍』って印が押してあったよ。そうか、マスが 有限なんで、溢れちゃうからお裾分けしますよって事だな。

コンピュータアーキテクチャ（オーム社）とOpenGLによる3Dグラフィック（朝倉書店）を頂いてきた。 こんな図書、今まで館内で見かけなかったぞ。今まで多分、蔵に入っていたんだろう。蔵が一杯に なったんで放出したんだな。町は貧乏なので、どこかの大学から横滑りしてきた本みたい。

両本とも教科書然としてました。演習問題がやけに多くて、ちょいと読みにくい。でも、無料なんで 許しちゃうぞ。もう15年も前の本ですけどね。アーキ本の付録に、CPUのクロック予測が出ていて、 それによると、2010年ぐらいは、10GHzぐらいになると予測されてました。

無邪気で実に宜しい。半導体を知らない人の予測だな。今じゃ、クロックスピードを上げて、CPU 性能を向上しようって魂胆はとっくに破綻してる。一人の天才よりも、よってたかって凡人に 仕事を割り振るのがいいって方針に転換してるものなあ。その最たるものとして、小さいやつは GPU、大きいやつは京なんだな。

代表的なワークステーションとして、今は亡きSunのウルトラ60が上げられていて、妙に懐かしく なっちゃいましたよ。

後は、演算器の話が面白かったり、アドレッシングで３アドレス、２アドレス、１アドレス方式が あるよとか。１アドレスの例としてスタックマシンが出てた。 詳細は、このあたりかな。世のVMはほとんどが、 １アドレス方式のスタックマシンだな。あれ？、複数のスタックを持ってるのは、1.5アドレス方式 になるんかな。

jacal

前回やった数式処理はPython流だったけど、やっぱり本場はLispだな。思い出したよ。微分の コードが出てたのは、 sicpなんだった。すっかり忘れていた。

;;;SECTION 2.3.2

(define (deriv exp var)
  (cond ((number? exp) 0)
        ((variable? exp)
         (if (same-variable? exp var) 1 0))
        ((sum? exp)
         (make-sum (deriv (addend exp) var)
                   (deriv (augend exp) var)))
        ((product? exp)
         (make-sum
           (make-product (multiplier exp)
                         (deriv (multiplicand exp) var))
           (make-product (deriv (multiplier exp) var)
                         (multiplicand exp))))
        (else
         (error "unknown expression type -- DERIV" exp))))
  :
;: (deriv '(+ x 3) 'x)
;: (deriv '(* x y) 'x)
;: (deriv '(* (* x y) (+ x 3)) 'x)

こんな方法が提示された。そして、結果の簡約化にも言及されたよ。以前はこの本を最後まで 読んでみよーって事で、理解が薄っぺらになってたけど、改めて読み直すと面白いな。

MITは、なんでこんな楽しい本を止めて、Pythonに蔵替えしちゃったんかな。何でも、教科書に 採用されたのは、Dive Into Pythonらしい。やっぱり、括弧に 抵抗を持たれたのかなあ。

Scheme 数式処理 でIEしてたら、Topにあの方が出てきた。 その後、進展はいかがでしょうか？ いきなり、Gaucheが下支えしてるページに飛んだけど、 最近、gaucheも新しい版が出たんで、更新しておこうかな。

そして、数値積分や微分を扱ったページも 見つけた。この方もSICPのファンのようです。遅延評価と言えば、Haskellだけど、Haskell方面で 面白そうなものはあるのかな？ 余りアカデミックな物より、使えるものの方が面白いのにね。 きっと、haskellは悟りを開くための修行用なんだな。

schemeで動く数式処理が無いか探して見つけたのが、The JACAL Symbolic Math System。 FreeBSDのportsになってたので、入れてみた。

[sakae@cdr ~]$ jacal
scm -ip1 -l /usr/local/lib/jacal/go
.............
JACAL version 1b7, Copyright 1989-2005 Aubrey Jaffer
JACAL comes with ABSOLUTELY NO WARRANTY; for details type `(terms)'.
This is free software, and you are welcome to redistribute it
under certain conditions; type `(terms)' for details.
;;; Type (math) to begin.
type qed; to return to scheme, type help; for help.
e0 :
;ERROR: "/usr/local/lib/jacal/sexp.scm": unbound variable:  tok:bump-column
; in expression: (tok:bump-column (string-length math:prompt) cip)
; in scope:
;   (echoing . #@define)
;   ()  procedure loop
;   (loop . #@let)
;   (cip obj . #@define)
;   (math:exit-cnt)  procedure <anon>
;   (math:prompt var-news-saved math:exit-saved . #@do)
;   ()  procedure batch1
; defined by load: "/usr/local/lib/jacal/sexp.scm"

;STACK TRACE
1; (#@define ((echoing (#@not (#@eq? (#@get-grammar (#@quote null ...
2; (#@define ((cip (#@current-input-port)) (obj #f))  (#@set! #@m ...
3; (#@do ((math:prompt #f math:prompt) (var-news-saved #@var-news ...
4; ((#@set-handlers!) (#@for-each (#@lambda (file) (batch (if (sy ...
5; (#@math)
6; (#@define ((hss (#@has-suffix? #@file (#@scheme-file-suffix))) ...
7; ((#@do-thunk (#@lambda () (#@cond (#@*syntax-rules* (require ( ...
8; ((#@case #@option #(#<unspecified> #f #\? #\: #\n #\u #\m #\s  ...
9; ((#@cond ((#@not #@*argv*) (#@set! #@*argv* (#@program-argumen ...

ははは、いきなり内臓を見せられた気分ですよ。こういう事もたまには有るのね。 何か足りないのかな。そんな風に思えるけど、debugの仕方も知らないし、あれこれ 悩むより、最新のやつを取ってくる方が簡単だろう。

取ってきて、全てのscmファイルを、/usr/local/lib/jacal内に押し込むと言う 暴挙に出てみた。

[sakae@cdr ~]$ jacal
scm -ip1 -l /usr/local/lib/jacal/go
.............
JACAL version 1c2, Copyright 1989-2005 Aubrey Jaffer
JACAL comes with ABSOLUTELY NO WARRANTY; for details type `(terms)'.
This is free software, and you are welcome to redistribute it
under certain conditions; type `(terms)' for details.
;;; Type (math) to begin.
type qed; to return to scheme, type help; for help.
e0 : diff(x^3 + 2 * x, x);

           2
e0: 2 + 3 x

ああ、動いた。後は、マニュアルに首ったけになって、どんな事が出来るか 試してみよう。案内によると、schemeとの間を行ったり来たり出来るみたい だし。面白そう。ああ、本当の面白さは、ソースの観賞なんですが、何処から 眼を付ければいいのだろうか。

maxima

Lisp方面で有名な数式処理システムがないかと探してみた。そしたら、 Maxima, a Computer Algebra Systemが出てきた。 FreeBSDではportsになってるので入れてみようとしたんだ。そしたら、豪華な事に、CLには 何を使いますか？だと。clisp、gcl、cmucl、sbclの中から選べですって。今回はたまたま入って いた、鉄鋼銀行コモンシスプを指定してみた。

そんでもって、使い方も探しておいたら、 Maxima による数式処理とか、PDFで出来た （大学の資料とおぼしきマニュアル）が沢山見つかったよ。一生懸命にmaximaを勉強させておいて、 試験には、maximaを使ってはいけませんとは、おかしな大学だなあ。どことは言わんけど。

早速起動してみる。

[sakae@cdr ~]$ maxima
Maxima 5.25.0 http://maxima.sourceforge.net
using Lisp SBCL 1.0.34
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) factor(2^32 - 1);
(%o1)                          3 5 17 257 65537
(%i2) factor(2^32 + 1);
(%o2)                             641 6700417
(%i3) factor(2**10 - 1);
(%o3)                               3 11 31
(%i4) factor(2**64 - 1);
(%o4)                    3 5 17 257 641 65537 6700417
(%i5) factor(2**127 - 1);
(%o5)               170141183460469231731687303715884105727

なかなか、メルセンヌ素数は見つからないですねぇ。上のマニュアルを見て行くと、グラフはgnuplotに お任せよーなんですね。新しいgnuplotだと、3Dグラフも使えるんだ。（昔やったのを思い出したぞ）

(%i6) to_lisp();

Type (to-maxima) to restart, ($quit) to quit Maxima.

MAXIMA> (* 3 11 31)

1023
MAXIMA> (to-maxima)
Returning to Maxima
(%o6)                                true
(%i7) :lisp (* 3 5 17 257 65537)

4294967295

lisp界隈へ行くのも勿論OK。手軽にやるなら、コロンコマンドも有ります。

(%i19) abc: 123;
(%o19)                                123
(%i20) :lisp $abc

123

maxima側の変数名は、lisp側では、頭に＄を付けて表される事は分かったけど、関数は、どう なっているんだろう？

(%i25) foo: a + b;
(%o25)                               b + a
(%i26) :lisp $foo

((MPLUS SIMP) $A $B)

こういうのは、Lisp側の片鱗が見えるんだけどな。やっぱりソースを見ないとだめかな？ その場合でも、こちらは 参考になるのかな？

(%i27) ? time;

 -- Function: time (%o1, %o2, %o3, ...)
     Returns a list of the times, in seconds, taken to compute the
     output lines `%o1', `%o2', `%o3', .... The time returned is
     Maxima's estimate of the internal computation time, not the
     elapsed time. `time' can only be applied to output line variables;
     for any other variables, `time' returns `unknown'.

     Set `showtime: true' to make Maxima print out the computation time
     and elapsed time with each output line.


  There are also some inexact matches for `time'.
  Try `?? time' to see them.

(%o27)                               true

あらかじめコマンド名が分かっていれば、オンラインヘルプが利用出来る。知らないコマンドは

(%i32) ?? prime;

 0: mod_big_prime  (Functions and Variables for zeilberger)
 1: next_prime  (Functions and Variables for Number Theory)
 2: prev_prime  (Functions and Variables for Number Theory)
 3: primep  (Functions and Variables for Number Theory)
 4: primep_number_of_tests  (Functions and Variables for Number Theory)
Enter space-separated numbers, `all' or `none':

上記のように、適当にあたりを付けて、？？ すれば出てくる。地理感のある人は、これだけで 十分だろう。しばし、いろんなキーワードを入れて、こんな事やあんな事まで出来ると感激に 浸ってしまったよ。

そんじゃ、timeを使って実行時間を計ってみる。題材は、Lispと言ったら階乗でしょう。

(%i37) 100000!$

(%i38) time(%o37);
(%o38)                              [5.196]

100000の階乗なんだけど、出力するととんでもない事になるんで、＄を使って出力を抑止。指示に 従って、実行時間を測ると、５秒でした。commonLispは、末尾再帰の最適化をやらないんで、こんな もんですかね？ そりゃ違うだろ。幾らなんでも、そんな深い再帰じゃスタックが足りなくて、 落ちちゃうはず。

そうすると、階乗をループで計算してるんだな。上記の時間のほとんどは、多倍長演算に費やされた 時間だろうな。

おまけ

数式をかっこよく表示させようと思って、Texmacsをいれようとした んだ。そしたら、LaTexのダウンロードが始まっちゃったんで速攻で中止した。もしコンパイルを 始めてしまったら、熱さでパソコンがダウンしちゃったかも。危ない、危ない！

で、代わりは無いものかと思って、ウブに入ってるOOを見たら、メニューには登録されていなかった。 普通の人には用が無いからなあ。コマンドラインから起動するんだな。

sakae@ubuntu:~$ openoffice.org -math

ドット疲れる、ドット付きなコマンドでしたよ。そんでもって、数学者が発明した記号は どうやって入力してくかだな。例によると、２次方程式の表示は、

x={-b+-sqrt{b^2-4 ac}}over {2 a}

らしい。下側の入力欄に上記の式を入力してくと、上段の表示欄に綺麗な式が描画されたよ。これには ちと感動しましただ。でも、こんな調子で数学者の発明した記号を入力するってのは、ちと苦痛。 もっと直感的な奴は無いかね？

探してみたら、Windows7限定らしいけど、 お絵かきが 出来るみたい。これって、ウィドウフォーン用の実験台だったんでしょうかね。

話は変わって、Ruby方面でも数式処理が無いかと探してみたら、 原先生のが見つかった。 もっと他にもあるかも知れないけど、親分が数学苦手って公言してたみたいで、数学方面は 流行らないのね。

Pythonだと、先に紹介した以外にも、 sageなんてのが見つかった。 そして、Pythonと統計ソフト R の間を取り持つ rsympyなんてのも見つかった。PythonはHubだなあ。