(define (test n ps is js) (define (up ps) (caddr (list-ref ps 5))) (define (left ps) (cadar ps)) (define (try *????) (for-each (lambda (x) (let ((p (car x)) (i (caadr x)) (j (cadadr x))) (if (not (member i is)) (test (+ n 1) (cons p ps) (cons i is) (cons j js))))) *????)) (case n ((0) (try *2??2)) ((1 2 3 4) (case (left ps) ((0) (try *2??0)) ((1) (try *2??1)) ((2) (try *2??2)))) ((5) (case (left ps) ((0) (try *22?0)) ((1) (try *22?1)) ((2) (try *22?2)))) ((6 12) (case (up ps) ((0) (try *0??2)) ((1) (try *1??2)) ((2) (try *2??2))))
((24) (begin (display ps) (newline) (display is) (newline) (display js)) (newline))))
<html> <head> <title>nth day after birth</title> <script language="JavaScript">
function quotient(a,b){
var r=a%b;
var q=(a-r)/b;
return q;}
var mon0=[0,31,59,90,120,151,181,212,243,273,304,334]; var mon1=[0,31,60,91,121,152,182,213,244,274,305,335];
function rd(y,m,d){
var a=(y-1)*365;
var b=quotient(y+3,4);
var c=quotient(y,400)-quotient(y,100);
var f=(((y%100)==0)?((y%400)==0):((y%4)==0))?mon1[m-1]:
mon0[m-1];
return a+b+c+f+d-1;}
function gr(rd){
var d0=rd+366-1;
var y400=quotient(d0,146097);
var d1=d0%146097;
var y100=quotient(d1/36524.25,1);
var d2=d1-146097+quotient(36524.25*(4-y100),1);
var y4=(y100==0)?quotient(d2,1461):
(24-quotient((36523-d2)/1460.96,1));
var d3=(y100==0)?(d2%1461):(d2-quotient(y4*1460.96,1));
var y1=((y100>0)&&(y4==0))?quotient(d3,365):
quotient(d3/365.25,1);
var d4=((y100>0)&&(y4==0))?d3%365:
(d3-1461+quotient(365.25*(4-y1),1));
var y=y400*400+y100*100+y4*4+y1;
var leap=((y%100)==0)?((y%400)==0):((y%4)==0);
var e=leap?mon1:mon0;
var mon=0;
for(i=0;e[i]<=d4;i++){mon=mon+1;}
var d=d4-e[mon-1]+1;
document.write(y.toString(10)+" "+mon.toString(10)+" "
+d.toString(10));
document.write("<br>");}
</script>
</head>
<body>
<script language="JavaScript">
document.write("<br>");
var s,i,y,m,d,n;
s = prompt("Type in your birthyear,birthmonth,birthday and n
as `2000,1,1,1000' without space", "");
i = s.indexOf(',');
y = s.substring(0,i);
s = s.substring(i+1,s.length);
i = s.indexOf(',');
m = s.substring(0,i);
s = s.substring(i+1,s.length);
i = s.indexOf(',');
d = s.substring(0,i);
n = s.substring(i+1,s.length);
document.write("y = " + y);document.write("<br>");
document.write("m = " + m);document.write("<br>");
document.write("d = " + d);document.write("<br>");
document.write("n = " + n);document.write("<br>");
y=y-0;
m=m-0;
d=d-0;
n=n-0;
var r=rd(y,m,d);
document.write("Your "+ n.toString(10) +"th day is ");
gr(r+n);document.write("<br>");
</script>
</body>
</html>
y = 1912 m = 6 d = 23 n = 15000 Your 15000th day is 1953 7 18
(define c (expt 10 100)) (do ((n 1 (+ n 2))) ((< (/ c n (expt 5 n)) 1) n)) => 141 (do ((n 1 (+ n 2))) ((< (/ c n (expt 239 n)) 1) n)) => 43これで必要な項数は判ったので, 105桁のbignumで計算する. aは tan-1(1/5), bはtan-1(1/239)で, 最後にaの16倍からbの4倍を引き, 最後の5桁を捨てる. opはループで毎回0と1に切り替わりそれに従って次の項を足したり引いたりする.
(define c (expt 10 105)) (define op 0) (define a 0) (do ((n 1 (+ n 2))) ((= n 143)) (set! a ((if (= op 0) + -) a (quotient (quotient c n) (expt 5 n)))) (set! op (- 1 op))) (define op 0) (define b 0) (do ((n 1 (+ n 2))) ((= n 45)) (set! b ((if (= op 0) + -) b (quotient (quotient c n) (expt 239 n)))) (set! op (- 1 op))) (quotient (- (* a 16) (* b 4)) 100000)で結果は
314159265358979323846264338327950288419716939937510 58209749445923078164062862089986280348253421170679πの100桁, 1000桁の値はすぐに見付かる. 私は城, 牧之内「計算機械」にあったのを知っているので書棚から取り出し同じであることを確認した.
852 0 11...110 負ストロブ
854 1 11...101
856 2 11...011
...
920 34 01...111
922 0 00...001 正ストロブ
924 1 00...010
926 2 00...100
...
990 34 10...000
9753←対応する奇数の素数
992 0 11...111 篩
994 1 11...111
996 2 11...111
...
1022 15 11...111
0番の篩の長語の各ビットはその上に書いてあるように右から3,5,7,9,...に対応する.従って最後の15番の篩の左端のビットは1121に対応する. EDSACのメモリー容量からすれば篩はもっと大きく出来るが, 篩の様子を水銀タンクで眺められるように, ちょうど1つのタンクの大きさにした.
篩は正のストロブを順に使い, 篩の右から順のビットとandをとって, 素数か合成数かをみる. 合成数なら次のストロブで次の篩のビットを調べに進む. 素数ならその素数をp とすると, 左の図のようにpビットごとにpの倍数があるから, 負のストロブを使っていま調べたビットからpビットおきに篩のビットを0に変えていく. 手順としてはこれだけだが, ストロブが下まで来たときに, 篩を次の長語にする; ストロブを先頭に戻す作業が必要である.
命令 オペランド アキュムレータ
47 A 34 @ C 992 D C 992 D
48 U 71 @ 71 C 992 D
49 S 843 F L F T 992 D
50 T 72 @ 72 T 992 D
51 A 33 @ H 922 D H 922 D
52 S 842 F P 70 F H 852 D
53 A 841 F P 6 F H 858 D
54 U 70 @ 70 H 998 D
55 S 840 F H 992 D P 6 F
56 G 69 @
47'で今素数を調べた篩を取り出し71'へ入れ, ファンクション部をTにした命令を72'へ入れる. 51'は素数を調べた正ストロブを乗数レジスタに入れる命令で, 篩を0にするには対応する負ストロブが必要である. それは7番地若い方にあるからまず70を引き, 素数のビット数だけ先から篩うので841番地にある素数の4倍の値を足し, 消すストロブの
番地を作る. 55'は負ストロブの範囲を超えたかのチェック. 69 T F 70 H 998 D 71 C 992 D 72 T 992 D 73 A 70 @ H 998 D H 998 D 74 G 53 @ 53 A 841 F P 6 F H1004 D 54 U 70 @ 70 H1004 D73'からは篩をクリアする命令を更新する. H 998 Dをもって53'へ. そこで再び素数の数だけストロボのアドレスを増やす. こうしている内に負ストロブの範囲を飛び出し57'へ進む. 58'で篩の番地を取り出し, アドレスを2増やし71'へ戻す. 72'のT命令も増やす. しかし篩の方が範囲を超える心配があるので, 63'でT1024 Dを引き, 篩の範囲が終わっていなければストロブの番地をもとへ戻す. それが66'からでストロブを取り出し, 番地から70を引いて1周戻し, 54'へいって次のサイクルに入る
e t a i n o s r l d h c u m f p y g w v b k x j q z 1 1 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 3 4 4 3 4 2孔が少ないとテープが丈夫なのか, ごみを減らそうとしたのか知らないが, Morseコードと同様な精神で出来ている. ITA2の特殊記号も欧文Morseコードと同じである.
P 0 01101 Q 1 11101 W 2 11001 E 3 10000 R 4 01010 T 5 00001 Y 6 10101 U 7 11100 I 8 01100 O 9 00011と数値に関係がない. そこで最上段の文字のコードを00000から01001に変更した. ただそうするとP 0のコードはまったく孔が開かず, ブランクテープと同じになってしまうので, テープでは一番左の位置はコードが0の時に孔が開くようにした.
1 ]=> (* (/ 99999 65536) .655364990234375) ;Value: .9999976144172251 1 ]=> (* (/ 1 65536) .655364990234375) ;Value: 1.00000761449337e-5だから, 純小数にして5桁出力する方針である.
この関係をグラフにすると, 上の図の上のようになる. 
8<1を読み込むと, 要素1には先行者8があったので, 1のcountを1増やし, 1を要素8のtopのリストにつなげる. したがってデータ構造はBのようになる.
キューを作るには, 先頭(Front)から末尾(Rear)へ至るリストを使う. キューの最後にはnilを入れておく. キューに要素を挿入するのは末尾の側で, 削除するのは先頭の側である. その場所を覚えるためにの変数fとrを用意し, キューが空の時は, f=r=nilとする(E). 
まずfの値0を出力, N=N-1とする. 0のtopからリストを辿り, それらに先行していたものがなくなった処理をする. つまり要素2のcountを1減らす. その結果countが0になったら, この要素も先行するものがなくなり, 出力できるようになったわけで, この要素をcountのキューにつなぐ(H). 0のcountは8のままだが, これはゴミである. Algorithm 2.2.3-Tでは, topのリストはそのままにしているが, 演習問題2.2.3-23にあるように, トポロジカルソートが出来なかったときの原因探索用に, 済んだリストはnilにしておくのがよい.
図Iは要素8も出力したところ, Jは2も出力したところである. キューの先頭をfで辿り出力を続ける. fがnilになったら終わりで, Nは0になったはずである. 0にならなければ, 元のデータがおかしかった(ループがあった)ので, topのリストを出力してみると, 原因がわかる.
(define (algorithm223t n ls)
(let* ((count (make-vector n 0)) (top (make-vector n '()))
(r '()) (f '()) (nn n))
(define (getrel ls) ;データを読み込む
(if (not (null? ls))
(let ((j (caar ls)) (k (cadar ls)))
(vector-set! count k (+ (vector-ref count k) 1))
(vector-set! top j (cons k (vector-ref top j)))
(getrel (cdr ls)))))
(define (outloop f) ;出力のループ
(define (eraserel p) ;topのリストを解消する
(if (not (null? p))
(let ((c (- (vector-ref count (car p)) 1)))
(if (> c 0) (vector-set! count (car p) c)
(begin (vector-set! count r (car p)) (set! r (car p))
(vector-set! count r '())))
(eraserel (cdr p)))
(begin (vector-set! top f '())
(outloop (vector-ref count f)))))
(if (not (null? f))
(begin (display f) (display " ") (set! nn (- nn 1))
(eraserel (vector-ref top f)))))
(getrel ls)
(do ((k 0 (+ k 1))) ((= k n)) ;count=0のキューを作る
(if (= (vector-ref count k) 0)
(begin (if (null? r) (set! f k) (vector-set! count r k))
(set! r k) (vector-set! count r '()))))
(outloop f) ;出力開始
(if (= nn 0) 'ok (list 'nn nn))))
(algorithm223t 9
'((8 1) (2 6) (6 4) (4 7) (7 5) (3 5) (0 2) (6 3) (8 4) (1 7)))
このプログラムを走らせると 0 8 2 1 6 3 4 7 5 が出力される.
A, B, C, などが記憶場所で, その上方にブロックを置き, そのブロックまでの距離を記憶している値とする. 従って書き込むにはその数値の分だけブロックを遠ざける手続きを繰り返す. また読み出すにはブロックを近づけ, 0になるまでの回数で知る.
つまり頂点0,1,2の面に対しては, それに隣る三角形の頂点5,10,3の三角形で, こういう面だけで出来る立体はそう容易には想像出来ない. 数学者の想像力にはいつもながら脱帽だ. かなりの日時をかけ, 苦心惨憺, 紆余曲折, といいつつも十分楽しんだ末, 出来た図はこれだ. 自画自賛すれば実にもっともらしい.
たしかに頂点5,10,3を結ぶ三角形は見える. 下は左が正二十面体の面, 右がそれに対応する大きい三角形である.
[2 1 0] [10 5 3] [3 2 0] [11 1 4] [4 3 0] [7 2 5] [5 4 0] [8 3 1] [1 5 0] [9 4 2] [2 10 1] [11 9 0] [11 10 2] [6 1 3] [3 11 2] [7 10 0] [7 11 3] [6 2 4] [4 7 3] [8 11 0] [8 7 4] [6 3 5] [5 8 4] [9 7 0] [9 8 5] [6 4 1] [1 9 5] [10 8 0] [10 9 1] [6 5 2] [7 8 6] [4 9 11] [8 9 6] [5 10 7] [9 10 6] [1 11 8] [10 11 6] [2 7 9] [11 7 6] [3 8 10]ためつすがめつ眺めているうちに, これは下の図のオレンジ色の正五角形の上に, 背の高い星形が乗っている; その正五角形が芯になる正十二面体であるらしいと分かって来た.
しかし頂点3に関わる部分だけで三角形は15もある. 都合180個の三角形の座標を計算しなければならぬ. 当然その計算の部分はSchemeで別にやることにした.
上が正二十面体である. その頂点0に隣る5個の頂点のなす正五角形は, 下の図にオレンジ色で示す面1,2,3,4,5になる. これを「0に隣る面」ということにしよう. (その辺には小さく0と付記した.)
さらに頂点2と3に隣る頂点の面を作ると
のように2に隣る面0,1,10,11,3と, 3に隣る面0,2,11,7,4が得られる,
ところで0に隣る面と2に隣る面は1と3で交わるからその交線は1から3へ引いた線である. 同様に0に隣る面と3に隣る面は2と4で交わるからその交線は2から4へ引いた線である. さらに2に隣る面と3に隣る面の交線は0から11へ至る.
従ってもとの正二十面体の面の1つの三角形0,2,3の中に3つの交線が現れる. その交点を下の図のようにPとする.
同様に三角形0,3,4の中に3つの交線が現れる. その交点を下の図のようにQとする.
この2点の座標が決まると, 大十二面体を描くことが出来る. それがこれだ. この正多面体はPやQの所が凹んでいて, 水が溜まる形である. 右下の頂点3を中心とした昔の陸軍の徽章のような星形がみられる. 3の星の1つの枝は0の星の枝でもあるから, この枝は2等辺三角形である.
ところでPやQの座標の計算法だが, 正二十面体をx軸方向から見た図を描くと
になる. この図で1, 10短い方の辺の長さを2とすると, 図の
0, 6の長さは1+√5. Oから0までは(1+√5)/2.
0, 3も1なので, これらからP, Qのz座標(3+√5)/2, (√5-1)/2が得られる.
これらを解くには (a+b√5)×(c+d√5)のような計算が頻発するから, この積がe+f√5であるとして, a,b,c,dからeとfを求めるプログラム
(define (mult a b c d) (list (+ (* a c) (* 5 b d)) (+ (* a d) (* b c))))を用意すると仕事がはかどった.
