還暦

　 Wikipediaによると

還暦（かんれき）とは干支（十干十二支）が一巡し誕生年の干支に還ること。

　だそうだ。
　つまり、12と10の最小公倍数ということか。

(lcm 12 10)
60

　うん。あってるみたいだ。

平面図が正方形で、正面図が円で、右側面図も円(3)

　首尾よく「平面図が正方形で、正面図が円で、右側面図も円」が描けたので気をよくした。
　よし。次は「平面図が正方形で、正面図が円で、右側面図は正方形」を描いてみようじゃないか。
　うーん……っていうか、それってただの円柱じゃね……orz。

平面図が正方形で、正面図が円で、右側面図も円(2)

　「はじめての数式処理ソフト―Maximaで楽しむ数式計算と物理グラフィック」のgnuplotのプログラムをつらつら読むと、どうやら正方形をZ軸に沿って積み上げているらしい。
　ということはまず、正方形を描けないといけないな。
　けれど単純に描くと、discreteを使ってこんな感じになってしまう。

plot2d([discrete,[[-1,-1],[-1,1],[1,1],[1,-1],[-1,-1]]],[x,-2,2],[y,-2,2],same_xy)$

　立体にするから三角関数で描く形にした方がよいのだけど

plot2d([discrete,[[cos(0),sin(0)],[cos(0.5*%pi),sin(0.5*%pi)],[cos(%pi),sin(%pi)],[cos(1.5*%pi),sin(1.5*%pi)],[cos(0),sin(0)]]],same_xy)$

　これでは、軌跡を描いているわけではないからだめだ。
　parametric指定で三角関数で正方形を描くには……うーん、こうか？

plot2d([[parametric,1,sin(t),[t,0,2*%pi]],
        [parametric,-1,sin(t),[t,0,2*%pi]],
        [parametric,cos(t),1,[t,0,2*%pi]],
        [parametric,cos(t),-1,[t,0,2*%pi]]],
       [x,-2,2],[y,-2,2],same_xy)$

　1辺を二回ずつ、描いているところがちょっと嫌だけど、まぁ、いいや。
　これを円柱のときのように立体化。

plot3d([[1,sin(t),h],[t,0,2*%pi],[h,-1,1],
        [-1,sin(t),h],[t,0,2*%pi],[h,-1,1],
        [cos(t),1,h],[t,0,2*%pi],[h,-1,1],
        [cos(t),-1,h],[t,0,2*%pi],[h,-1,1]],
       [x,-2,2],[y,-2,2],[z,-2,2],same_xyz)$

　Z軸の位置によってかわる正方形の一辺の長さは三平方の定理からわかるから

半径 = 1
長さ = 2*sqrt(1-z^2)

　その分、X,Y座標を調整してやれば、たぶん……。

plot3d([[sqrt(1-h^2),sin(t)*sqrt(1-h^2),h],[t,0,2*%pi],[h,-1,1],
        [-sqrt(1-h^2),sin(t)*sqrt(1-h^2),h],[t,0,2*%pi],[h,-1,1],
        [cos(t)*sqrt(1-h^2),sqrt(1-h^2),h],[t,0,2*%pi],[h,-1,1],
        [cos(t)*sqrt(1-h^2),-sqrt(1-h^2),h],[t,0,2*%pi],[h,-1,1]],same_xyz)$

　おっ、できた。

　平面図

　正面図

　右側面図

　すばらしい。

　そして、45度からの側面

　実はこの角度から見たとき、楕円じゃないか、と思ってたしかめたかったのだ。
　きっとこれは長軸、正方形の対角線、短軸、正方形の一辺(=円の直径)の楕円だよね。ほー。

平面図が正方形で、正面図が円で、右側面図も円(1)

　竹内薫「はじめての数式処理ソフト―Maximaで楽しむ数式計算と物理グラフィック」の中の、平面図が正方形で、正面図が円で、右側面図も円の立体に関する問題。
　そこに「Maximaでは複数の3次元グラフを重ねられない」とある。
　——ほんと？
　と思ったのも、本の中ではgnuplotで図を描いているのだけど、それって変じゃね？ と思ったからだ。Maximaもgnuplotで描いているんだから描けそうな気がする。どうなんだろ？

　「平面図が正方形で、正面図が円で、右側面図も円」ということは横から見ると、正方形の円柱を寝転がして二個、90度ちがいで重ねれば、よいような……でもそもそも円柱ってどうやれば、いいんだろ？
　まずは円か。
　これはいちおう、描ける。

plot2d([parametric,cos(t),sin(t),[t,0,2*%pi]],same_xy)$

　で、これに高さを与えれば、いいんだから……こうですかね？

plot3d([cos(t),sin(t),h],[t,0,2*%pi],[h,-1,1],same_xy)$

　で、これを転がす、と。

plot3d([h,cos(t),sin(t)],[t,0,2*%pi],[h,-1,1],same_xy)$

　さらにこれを90度、ずらして二個、描くと。

plot3d([[h,cos(t),sin(t)],[t,0,2*%pi],[h,-1,1],[cos(t),h,sin(t)],[t,0,2*%pi],[h,-1,1]],same_xy)$

　——えっ？
　やっぱ、3次元グラフを重ねて描けるじゃないか。
　まぁ、10年前の本だから当時は描けなかったのかもしれないけれど。

　平面図

　正面図
　——このはみ出ているところを削ることができれば、答えなんだけどなぁ。これはどうすれば、いいんだろう？ まぁ、とりあえず、3次元グラフを重ねられるのなら「はじめての数式処理ソフト―Maximaで楽しむ数式計算と物理グラフィック」のgnuplotのプログラムもmaximaできっと書けるんじゃね？
　こうかな？

f(z,t):=t*sqrt(1.0-z^2);

plot3d([[f(-0.9,u),f(-0.9,v),-0.9],[u,-1,1],[v,-1,1],
        [f(-0.8,u),f(-0.8,v),-0.8],[u,-1,1],[v,-1,1],
        [f(-0.7,u),f(-0.7,v),-0.7],[u,-1,1],[v,-1,1],
        [f(-0.6,u),f(-0.6,v),-0.6],[u,-1,1],[v,-1,1],
        [f(-0.5,u),f(-0.5,v),-0.5],[u,-1,1],[v,-1,1],
        [f(-0.4,u),f(-0.4,v),-0.4],[u,-1,1],[v,-1,1],
        [f(-0.3,u),f(-0.3,v),-0.3],[u,-1,1],[v,-1,1],
        [f(-0.2,u),f(-0.2,v),-0.2],[u,-1,1],[v,-1,1],
        [f(-0.1,u),f(-0.1,v),-0.1],[u,-1,1],[v,-1,1],
        [f(-0.0,u),f(-0.0,v),-0.0],[u,-1,1],[v,-1,1],
        [f(0.9,u),f(0.9,v),0.9],[u,-1,1],[v,-1,1],
        [f(0.8,u),f(0.8,v),0.8],[u,-1,1],[v,-1,1],
        [f(0.7,u),f(0.7,v),0.7],[u,-1,1],[v,-1,1],
        [f(0.6,u),f(0.6,v),0.6],[u,-1,1],[v,-1,1],
        [f(0.5,u),f(0.5,v),0.5],[u,-1,1],[v,-1,1],
        [f(0.4,u),f(0.4,v),0.4],[u,-1,1],[v,-1,1],
        [f(0.3,u),f(0.3,v),0.3],[u,-1,1],[v,-1,1],
        [f(0.2,u),f(0.2,v),0.2],[u,-1,1],[v,-1,1],
        [f(0.1,u),f(0.1,v),0.1],[u,-1,1],[v,-1,1],
        [f(0.0,u),f(0.0,v),0.0],[u,-1,1],[v,-1,1],
        [f(u,1),f(u,1),u],[u,-1,1],[v,-1,1],
        [f(u,1),f(u,-1),u],[u,-1,1],[v,-1,1],
        [f(v,-1),f(v,1),v],[u,-1,1],[v,-1,1],
        [f(v,-1),f(v,-1),v],[u,-1,1],[v,-1,1]])$

　あかん、むっちゃ、重い。

円

　今だにGrapherではできた

1 > abs(x - y)

　のグラフは描けないけれど、そもそも円はどうやって描くのだろう。
　Grapherなら

　であっさり描けた。円。
　ところがgnuplot的にも簡単ではないのか、三角関数を使う方法しか、見当たらない。ちなみに

plot2d(sqrt(2^2-x^2),[x,-5,5],[y,-5,5],same_xy)$

　では半円になってしまった。

　いろいろ調べてみたらparametricを使えば、いいらしい。三角関数を使うやり方だけど。
　で、最終的にinfoに書かれていたのを発見して脱力する。それを元に描く。

plot2d([[parametric, cos(t), sin(t), [t,0,2*%pi]]],same_xy)$

　ここからおれが得るべき教訓は、まずinfoを見ろ、だな。

　三角関数を使わないのなら……半円にはなるわけだから……。
　solve()してみる。

(%i56) solve(1=sqrt(x²+y²),y);

　ということは

plot2d([sqrt(2^2-x^2),-sqrt(2^2-x^2)],[x,-5,5],[y,-5,5],same_xy)$

　で、どうよ。

　まぁ、こんなもんか。

　と思っていたのだけど、implicit_plotというパッケージを使えば、いちおう描けるみたいだ。

(%i102) load(implicit_plot)$

(%i103) implicit_plot(sqrt(x^2 + y^2) = 2, [x, -5, 5], [y, -5, 5],same_xy)$

(%i104) 

imaxima

　FreeBSDでimaximaが動いているせいでアスキーアートのmaximaに我慢できなくなった。
　MacBook AirにTeXをいれることにする。
　今度は用心深く最低限のTeX。

brew install basictex

　latexが見えないというのでパスを設定、そのほか。

;; imaxima
(setq exec-path (cons "/usr/local/texlive/2019basic/bin/x86_64-darwin" exec-path))
(autoload 'imaxima "imaxima" "Image support for Maxima." t)
(autoload 'imath-mode "imath" "Interactive Math minor mode." t)

(setq imaxima-fnt-size "large")
(setq imaxima-use-maxima-mode-flag t)

　やっと数式がきれいになった。グラフのインライン表示もできるようになった。

　ただ、Grapherではできた

1 > abs(x - y)

　のグラフはどうやったら描けるんだろう？ できないのかな？

Maxima on Emacs

　同梱されているmaxima.elだと、数式はアスキーアートになる。

(%i3) expand ((x + y)^6);
           6        5       2  4       3  3       4  2      5      6
(%o3)     y  + 6 x y  + 15 x  y  + 20 x  y  + 15 x  y  + 6 x  y + x

　それをTeXで

　という風に、きれいに表示するimaxima.elもあるので、TeXをインストールした。そうしたらひどい目にあった。
　MacTeXをいれたのだが、ディスクが満杯になってしまったのだ。MacBook Airの動きがおかしくなる。いや、もしかしたらMacTeXのせいではなく、タイミングわるく別の何かがあったのかもしれないのだが——。
　結局、MacTeXインストール以前の過去にもどることにした。
　TimeCapsuleを使って。
　——24時間以上かかった。24時間前にもどるのに24時間以上かかるとは。
　はあ。

　wxMaximaではなぜか、gnuplotがうまく動かないし。突然、ロシア語になるし。
　しかたないのでアスキーアートで使うか、といじっていたらdescribe()が効かないことに気づいた。

(%i33) describe(sum);
No exact match found for topic `sum'.
Try `?? sum' (inexact match) instead.

(%o33) false

　——はい？
　試しにFreeBSDへいれてみたMaximaはちゃんと動く。imaximaもちゃんと動く。数式もきれいだ。もうFreeBSDでいいじゃないか、と思ったけれど、CPUパワーがちがいすぎるし、ネットごしだし、オフラインだと使えない。
　SBCLが見えるところにないことが原因か、と

brew install sbcl

　したり、

brew install --build-bottle maxima

　buildしなおしたりしたけれど、解決せず。
　FreeBSDでは出力されていないので原因は起動直後にでるこのwarningらしいのだが。

warning: ignoring an empty documentation index in /usr/local/Cellar/maxima/5.43.0/share/info/./

　わからん。なぜだ。
　Maximaでswankを動かしてslimeで裏から入って調べてみたりして——ソースを見てみるか、とbuildしたときにもってきたソースをのぞいたらinfoディレクトリの下にあるmaxima-index.lispのサイズが妙にちがう。いかにもな名前だし、このlispが動くことによってdescribe()のデータベースができるのだ、と当たりはつけていたのだけど。
　中身をコピーしてみた。
　立ち上げなおした。
　——動いた。

(%i1) describe(sum);


 -- Function: sum (<expr>, <i>, <i_0>, <i_1>)

     Represents a summation of the values of <expr> as the index <i>
     varies from <i_0> to <i_1>.  The noun form ''sum' is displayed as
     an uppercase letter sigma.

【中略】

(%o1) true

　くそお。

Maxima

　Maximaを入れた。

brew install maxima

　中身はCommonLispで書かれている、と聞いて興味を引かれたからだ¹。元々、Mathematicaがバンドルされていることを知ってRaspberry Piを買うくらいには数学に憧れている文系だ。当然、インストールする。それがCommonLispでできているとなればなおのこと。数学は全然、わからないけれど。

Maxima 5.43.0 http://maxima.sourceforge.net
using Lisp SBCL 1.5.3
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.

　installしていてあれ？ と思った。
　SBCLがインストールされなかった²。必要なんじゃねーの？

brew info maxima

　してみると、

==> Dependencies
Build: sbcl ✘
Required: gettext ✔, gnuplot ✔, rlwrap ✔

　となっている。あ、requireではなく、SBCLでbuildされているのか。本体はコアダンプだから必要ないというわけか。

Footnotes:

¹

Emacsから使えるということも一因。

²

roswellには入っているけれど、brew的には入っていない。

Grapher

　奥の方にしまいこまれるようにMacBook Airには「Grapher」というアプリケーションがある。数式を指定するとグラフを描いてくれる。

y=x

　とか。

　じゃ、こんなんは？

1 > abs(x - y)

　描けた。しかも範囲も認識している。
　当然、3Dも描ける。Exampleにはこんなんまである。

　しかもマウスでぐりぐり動く。
　すげー。数式の意味はさっぱりだけど。

必勝法

>> モンテカルロシミュレーション

　当該ページで倍増法を必勝法みたいにいっているけど、ほんとうなのかなぁ。
　単純に考えて

1. ゲームのルールを決めることができるハウスがゲストの必勝法(ハウスにとって必敗法)があるゲームを採用するはずがない。
   
2. クラップスというゲーム自体、長いこと、存在する。
   
3. ということは必勝法が存在しない、すくなくとも今だ、見つけられていない。
   

　ということになるんじゃないか。

　いちおう自分でもモンテカルロしてみる。

　1000回試行

　あれ？
　なんか、いけそうに見える。
　540回ぐらいの大きな負けを超えることができれば、あとは順調に勝ちが積み上がっていくような……。

　10000回試行

　100000回試行

　1000000回試行

　10000000回試行

　——なんなんだ、これは。フラクタルかよ！
　どこまでいっても大負けになる瞬間がある。10000000回試行では瞬間的に1200000000円——12億……。どんなに勝ちを積んでもフッ飛んでしまうことがあるのか。おそろしい。

クラップス

　まず、クラップスのルールを説明します。クラップスは２個のサイコロを用いておこないます。２個の目の合計は２から１２までの値をとりま す。（３．５節で示したとおり、７の出る確率が最も大きい。）は２個のサイコロを振って、その和が７または、１１の時（これをナチュラルといいま す）は勝ちで,２は１２の時（これをクラップスといいます）は負けです。もしそれ以外のとき、すなわち、４、５、６、は、１０がでた ときは、その値をポイントといい、ひきつづき同じポイントまたは７が出るまで、サイコます。そして同じポイントがでれば、プレーヤの勝ち、それ以 前に７が出れば、負けとなります。もしプレーヤが勝てば、自分の賭け金の２倍の払い戻しを受けます。
　このゲームは、たいへんよくできて、ごく僅かな確率の差で、胴元が勝つようにできています。

モンテカルロシミュレーション

　確率を計算してみる。
　まず、合計値の分布は次のようになる。

合計値	回数
2	1
3	2
4	3
5	4
6	5
7	6
8	5
9	4
10	3
11	2
12	1

　次に合計値ごとの確率を求めてみる。

合計値	回数	勝ち	負け
2	1	0	1 / 36
3	2	0	2 / 36
4	3	(3 / 36) * (3 / (3 + 6))	(3 / 36) * (6 / (3 + 6))
5	4	(4 / 36) * (4 / (4 + 6))	(4 / 36) * (6 / (4 + 6))
6	5	(5 / 36) * (5 / (5 + 6))	(5 / 36) * (6 / (5 + 6))
7	6	6 / 36	0
8	5	(5 / 36) * (5 / (5 + 6))	(5 / 36) * (6 / (5 + 6))
9	4	(4 / 36) * (4 / (4 + 6))	(4 / 36) * (6 / (4 + 6))
10	3	(3 / 36) * (3 / (3 + 6))	(3 / 36) * (6 / (3 + 6))
11	2	2 / 36	0
12	1	0	1 / 36
総計	36	244/495	251/495

　検算：244/495 + 251/495 => 1

　これは毎回1円ずつ、495回、賭けたら7円のマイナスになるということだから、回収率は(495 - 7)/495*100で、約98.58585858585859％ということになる。

　競馬よりもむちゃくちゃましだけど、いずれにしても負ける。

ユークリッドの互除法

　一松信の「暗号の数理」を読んでいたら「ユークリッドの互除法」というものがでてきた。二つの正の整数m、nの最大公約数を求めるアルゴリズムである。下記のような手順だ。

1. mとnを比較し、必要なら入れかえてm≧nとせよ。
   
2. mをnで割り、商をq、余りをrとせよ。
   
3. 余りrが0(割り切れた)なら、そのときの除数nが最大公約数である。これで完了。
   
4. 余りrが0でなければ、除数nと余りrとをm、n,と置き直して、2.へ戻れ、
   

   この操作を反復すれば、いつか必ず3.に逹っして、完了するであろう。
   

   　——一松信「暗号の数理」より
   

　で、なんで、これで最大公約数が求まるのだろう……。
　ピンとこなくて考えこんでしまった。
　うーむ。
　n自身が最大公約数のとき、求まるのはわかる。では2周目はどうだろう。一回目のrが最大公約数の場合。もちろんそのときはrでnが割り切れるはずだ。もちろん、mも割り切れる(mはnとrに分解できるのだから)。
　あ、なるほど。そういうことか。
　あとは以下同文。

(defun Euclid(m n)
  "ユークリッドの互除法。m、nの最大公約数を求める。"
  (if (< m n)
      (Euclid n m)
    (progn
      (let ((q (/ m n))
     (r (mod m n)))
 (if (= r 0)
     n
   (Euclid n r))))))

(Euclid 25 5)
5
(Euclid 7 5)
1