十進BASIC 第2掲示板過去ログ2010

LET t0=TIME

LET N=24 !位数

LET w=INT((N-1)/2)+1 !概算
DIM A(w) !１～ｎまでの奇数列、偶数列

LET m=0
FOR i=1 TO N
   IF MOD(i,2)=1 THEN !奇数なら
   !IF MOD(i,2)=0 THEN !偶数なら　←←←←←
      LET m=m+1
      LET A(m)=i
   END IF
NEXT i
redim A(m)


LET ANSWER_COUNT=0

DIM B(m)
FOR i=0 TO 2^(m-1)-1 !ビットパターンで検証する　※MSB=0
   MAT B=CON
   LET t$=right$(REPEAT$("0",m)&BSTR$(i,2),m) !m桁
   FOR k=1 TO LEN(t$)
      IF t$(k:k)="1" THEN LET B(k)=-1 !0:1、1:-1へ
   NEXT k

   IF DOT(A,B)=0 THEN !±1*1 + ±1*3 + … + ±1*(2*m+1)を計算する
      LET ANSWER_COUNT=ANSWER_COUNT+1
      !!!PRINT "DATA ";CHR$(34);t$;CHR$(34)
   END IF
NEXT i

PRINT "個数=";ANSWER_COUNT !奇数の個数
PRINT


PRINT "計算時間=";TIME-t0

END

100 LET N=100
110 LET c=0
120 FOR i=2 TO N
130    FOR k=2 TO i-1 !約数を確認する
140       IF MOD(i,k)=0 THEN GOTO 170 !ひとつでも割り切れるなら、素数でない
150    NEXT k
160    LET c=c+1 !素数
170 NEXT i
180 PRINT c !結果を表示する
190 END

SET TEXT font "Courier New Bold Italic",24
PLOT TEXT ,AT 0.3,0.2: "ABCabc"

PLOT TEXT ,AT 0.1,0.8: "ABCabcあいう漢字" !元の書体と大きさ

DRAW PlotText("ＭＳ 明朝","",24,"ABCabcあいう漢字") WITH SHIFT(0.1,0.7)
DRAW PlotText("ＭＳ 明朝","太字",0,"ABCabcあいう漢字") WITH SHIFT(0.1,0.6)
DRAW PlotText("","斜体",0,"ABCabcあいう漢字") WITH SHIFT(0.1,0.5)
DRAW PlotText("ＭＳ 明朝","太字 斜体",24,"ABCabcあいう漢字") WITH SHIFT(0.1,0.4)

END

EXTERNAL PICTURE PlotText(f$,a$,s,s$) !原点を基準に、太字（Bold）や斜体（Italic）で文字を描く
IF f$<>"" OR s<>0 THEN SET TEXT font f$,s

IF POS(a$,"斜体")>0 THEN LET a=0.5 ELSE LET a=0
DRAW TEXT(s$) WITH SHEAR(a)

LET dx=worldx(pixelx(0)+1) !１ドットずらす
LET dy=worldy(pixely(0)+1)
IF POS(a$,"太字")>0 THEN DRAW TEXT(s$) WITH SHEAR(a)*SHIFT(dx,0)
!IF POS(a$,"太字")>0 THEN DRAW TEXT(s$) WITH SHEAR(a)*SHIFT(0,dy) !※必要に応じて
!IF POS(a$,"太字")>0 THEN DRAW TEXT(s$) WITH SHEAR(a)*SHIFT(dx,dy) !※必要に応じて
END PICTURE

EXTERNAL PICTURE TEXT(s$) !原点を基準に文字を描く
PLOT TEXT ,AT 0,0: s$
END PICTURE

!平面上のベクトル方程式とそのグラフ

OPTION ARITHMETIC COMPLEX

DEF v(a,b)=COMPLEX(a,b) !複素数の和差、実数倍の計算を対応させる

DEF fnDOT(a,b)=( a*conj(b) + conj(a)*b ) / 2 !内積　a1*b1+a2*b2
!絶対値　ベクトル |a|^2=a・a　　複素数 |z|^2=z*conj(z)


SET WINDOW -8,8,-8,8 !表示領域を設定する
DRAW grid !XY座標
ASK PIXEL SIZE (-8,-8; 8,8) w,h !画像の縦横の大きさ(ドット単位)を調べる

LET cEps=0.1 !精度　※調整が必要である

SET POINT STYLE 1 !点の形状


!例　点A(1,0)を通り、方向ベクトル(2,3)である直線の方程式を求めよ。
!参考　複素数による表現
!　点A（α=a+b*i）、d（β=c+d*i）とすると
!　conj(β)*z-β*conj(z)=conj(β)*α-β*conj(α)

LET OA=v(1,0)
LET d=v(2,3)

FOR t=-3 TO 3 !t=[-∞,∞]　※範囲は調整が必要である
   LET OP=OA+t*d !(x,y)=(1,0)+t*(2,3)=(1+2*t,3*t)
   PLOT LINES: Re(OP),Im(OP);
NEXT t
PLOT LINES


!例　２点A(2,5)、B(-1,3)を結ぶ線分ABの方程式を求めよ。
!参考　複素数による表現
!　点A（α=a+b*i）、点B（β=c+d*i）とすると、直線ABは
!　(conj(β)-conj(α))*z+(β-α)*conj(z)=conj(β)*α-β*conj(α)

LET OA=v(2,5)
LET OB=v(-1,3)

FOR t=0 TO 1 !t=[0,1]
   LET OP=(1-t)*OA+t*OB !(x,y)=(1-t)*(2,5)+t*(-1,3)
   PLOT LINES: Re(OP),Im(OP);
NEXT t
PLOT LINES


!例　点A(2,5)、ベクトルn=(4,3)に垂直な直線の方程式を求めよ。（内積を用いた表現）
!参考　複素数による表現
!　点A（α=a+b*i）、n（β=c+d*i）とすると
!　conj(β)*z+β*conj(z)=conj(β)*α+β*conj(α)

LET a=v(2,5)
LET n=v(4,3)

FOR j=0 TO h !画面全体を走査する
   LET y=worldy(j) !ドットをxy座標に変換する
   FOR i=0 TO w
      LET x=worldx(i)

      LET p=v(x,y) !n・(p-a)=0となる点Pの軌跡
      IF ABS( fnDOT(n,p-a) )<cEps THEN PLOT POINTS: x,y

   NEXT i
NEXT j


!例　点C(-1,1)を中心とする半径3の円の方程式を求めよ。
!参考　複素数による表現
!　点C（α=a+b*i）とすると、|z-α|=r

LET OC=v(-1,1)
LET r=3

FOR j=0 TO h !画面全体を走査する
   LET y=worldy(j) !ドットをxy座標に変換する
   FOR i=0 TO w
      LET x=worldx(i)

      LET OP=v(x,y) !|p-c|=rとなる点Pの軌跡
      IF ABS( ABS(OP-OC)-r )<cEps THEN PLOT POINTS: x,y

   NEXT i
NEXT j



!例　２点A(-5,-3)、B(2,4)を直径とする円の方程式を求めよ。（内積を用いた表現）
!参考　複素数による表現
!　点A（α=a+b*i）、点B（β=c+d*i）とすると
!　arg((z-α)/(z-β))=±PI/2　または　|z-(α+β)/2|=|α-β|/2

LET OA=v(-5,-3)
LET OB=v(2,4)

FOR j=0 TO h !画面全体を走査する
   LET y=worldy(j) !ドットをxy座標に変換する
   FOR i=0 TO w
      LET x=worldx(i)

      LET OP=v(x,y) !PA・PB=0（円周角の定理）となる点Pの軌跡
      IF ABS( fnDOT(OA-OP,OB-OP) )<cEps THEN PLOT POINTS: x,y

   NEXT i
NEXT j


!例　点O(0,0)、A(2,0)、B(1,2)として、点PをOP=s*OA+t*OB、実数s,tとする。
!　①1≦s+t≦3、0≦s、0≦tの範囲を動くとき、点Pの存在範囲を図示せよ。
!　②2*s+t=1を満たすとき、点Pの存在範囲を図示せよ。

LET OA=v(2,0)
LET OB=v(1,2)

FOR s=0 TO 3 STEP 2^(-6) !s≦3　※刻みは調整が必要である
   FOR t=0 TO 3 STEP 2^(-6) !0≦sより、s+t≦3≦3+sとなる。これより、t≦3
      IF 1<=s+t AND s+t<=3 THEN !①
         LET OP=s*OA+t*OB
         PLOT POINTS: Re(OP),Im(OP)
      END IF
   NEXT t
NEXT s

FOR s=-4 TO 4 STEP 2^(-6) !t=[-∞,∞]　※範囲と刻みは調整が必要である
   LET t=1-2*s !2*s+t=1より
   LET OP=s*OA+t*OB
   PLOT POINTS: Re(OP),Im(OP) !②
NEXT s


!例　Ｘ軸上を転がる半径1の円の円周上の点Pの軌跡（サイクロイド）

LET r=1
FOR t=-3*PI TO 3*PI STEP 0.1 !※範囲と刻みは調整が必要である
   LET OA=v(r*t,r) !中心A、Ｘ軸との接点Bとすると、OB=弧BPより
   LET AP=v(-r*SIN(t),-r*COS(t))
   LET OP=OA+AP !点Pの軌跡
   PLOT LINES: Re(OP),Im(OP);
NEXT t
PLOT LINES


END

!直列多重振り子
!参考サイト　http://www.aihara.co.jp/~taiji/pendula-equations/present-node5.html

LET N=3 !振り子の数

DIM m(N) !振り子のおもりの質量
DATA 0.2, 0.3, 0.2
MAT READ m

DIM l(N) !振り子の糸の長さ
DATA 3, 3, 4
MAT READ l

DIM th(N) !振り子のy軸に対する角度
LET th(1)=2*PI/3 !初期値
LET th(2)=-PI/2
LET th(3)=PI/6

!θ[d]（d=1,n）に関するラグランジュ運動方程式
!   m[d,n]*l[d]^2*d/dt^2{θ[d]}
! + Σ[i=1,d-1] m[d,n]*l[d]*l[i]*d/dt^2{θ[i]}*C[d,i]
! + Σ[i=d+1,n] m[i,n]*l[d]*l[i]*d/dt^2{θ[i]}*C[d,i]
! + Σ[i=1,d-1] m[d,n]*l[d]*l[i]*(d/dt{θ[i])^2*S[d,i]
! + Σ[i=d+1,n] m[i,n]*l[d]*l[i]*(d/dt{θ[i])^2*S[d,i]
! + m[d,n]*g*l[d]*SIN(θ[d])
! = 0

LET G=9.8 !重力加速度

FUNCTION sm(k,l) !m[k,l]=Σ[j=k,l] mj
   local s,j
   LET s=0
   FOR j=k TO l
      LET s=s+m(j)
   NEXT j
   LET sm=s
END FUNCTION
DEF C(i,j)=COS(th(i)-th(j)) !C[i,j]=COS(θi-θj)
DEF S(i,j)=SIN(th(i)-th(j)) !S[i,j]=SIN(θi-θj)

!２階微分方程式を連立１階微分方程式にする
!　ωd=d/dt{θd}とすると、R*d/dt^2{θd}=vよりd/dt{ωd}=[INV(R)*v]d

DIM R(N,N),v(N)
SUB FNC(t,x(), f()) !連立微分方程式 d/dt{Xi}=f(t,Xi)
   FOR d=1 TO N !行
      FOR i=1 TO N !列
         IF i>d THEN !右上部分
            LET R(d,i)=sm(i,N)*l(d)*l(i)*C(d,i)
         ELSEIF i=d THEN !対角線
            LET R(d,i)=sm(d,N)*l(d)^2
         ELSE !左下部分
            LET R(d,i)=sm(d,N)*l(d)*l(i)*C(d,i)
         END IF
      NEXT i
   NEXT d
   !!!MAT PRINT R;

   FOR d=1 TO N !行
      LET ss=0
      FOR i=1 TO d-1
         LET ss=ss-sm(d,N)*l(d)*l(i)*x(i)^2*S(d,i)
      NEXT i
      FOR i=d+1 TO N
         LET ss=ss-sm(i,N)*l(d)*l(i)*x(i)^2*S(d,i)
      NEXT i
      LET ss=ss-sm(d,N)*G*l(d)*SIN(th(d))

      LET v(d)=ss
   NEXT d
   !!!MAT PRINT v;

   MAT R=INV(R)
   MAT f=R*v
END SUB


SET WINDOW -10,10,-12,8 !描画領域

DIM w(N) !振り子のy軸に対する角速度
MAT w=ZER

LET h=0.01 !時間刻み幅 ⊿t
FOR t=0 TO 30 STEP h !※調整すること

   SET DRAW mode hidden !ちらつき防止開始
   CLEAR
   DRAW PendulumN(1) !親から順に描画する
   SET DRAW mode explicit !ちらつき防止終了

   !WAIT DELAY 0.2 !※必要に応じて

   !４次のルンゲ・クッタ法で連立微分方程式 d/dt{Xi}=f(t,Xi)を解く
   DIM k1(N),k2(N),k3(N),k4(N), f(N)

   DIM x(N),TT(N)
   MAT x=w
   CALL FNC(t,x, f)
   MAT k1=h*f

   MAT TT=(1/2)*k1
   MAT x=w+TT
   CALL FNC(t+h/2,x, f)
   MAT k2=h*f

   MAT TT=(1/2)*k2
   MAT x=w+TT
   CALL FNC(t+h/2,x, f)
   MAT k3=h*f

   MAT x=w+k3
   CALL FNC(t+h,x, f)
   MAT k4=h*f

   FOR i=1 TO N
      LET w(i)=w(i)+(k1(i)+2*(k2(i)+k3(i))+k4(i))/6 !x(t+h)=x(t)+(k1+2*k2+2*k3+k4)/6
   NEXT i

   MAT TT=h*w !ωd=d/dt{θd}より
   MAT th=th+TT
   !!!MAT PRINT th;

NEXT t


PICTURE PendulumN(i) !ｎ重の振り子を描く
   LET LL=l(i)
   LET AA=th(i)
   DRAW Pendulum(LL,m(i)) WITH ROTATE(AA) !親
   IF i<N THEN DRAW PendulumN(i+1) WITH SHIFT(LL*SIN(AA),-LL*COS(AA)) !階層関係（子へ）
END PICTURE

PICTURE Pendulum(L,r) !原点を基準に振り子を描く
   PLOT LINES: 0,0; 0,-L !糸
   DRAW disk WITH SCALE(r)*SHIFT(0,-L) !おもり
END PICTURE

END

!図形とベクトル方程式（三角形の心）

!平面の点を「平面のベクトルとみる」と「複素数とみる」と考えられる

OPTION ARITHMETIC COMPLEX

DEF v(a,b)=COMPLEX(a,b) !複素数の和差、実数倍の計算を対応させる

DEF fnDOT(a,b)=( a*conj(b) + conj(a)*b ) / 2 !内積　a1*b1+a2*b2
!絶対値　ベクトル |a|^2=a・a　　複素数 |z|^2=z*conj(z)

DEF fnNormalize(a)=a/ABS(a) !正規化

!------------------------------ ここまでがサブルーチン


SET WINDOW -8,8,-8,8 !表示領域を設定する
DRAW grid !XY座標
ASK PIXEL SIZE (-8,-8; 8,8) w,h !画像の縦横の大きさ(ドット単位)を調べる

LET cEps=0.1 !精度　※調整が必要である

SET POINT STYLE 1 !点の形状


!例　三角形ABCの心
!　　　C
!　b ／＼ a
!　A ── B
!　　 c

LET OA=v(-4,-5)
LET OB=v(3,-4)
LET OC=v(4,6)

LET AB=OB-OA !辺AB
LET BC=OC-OB !辺BC
LET CA=OA-OC !辺CA

PLOT LINES: Re(OA),Im(OA); Re(OB),Im(OB) !辺ABを描く
PLOT LINES: Re(OB),Im(OB); Re(OC),Im(OC) !辺BC
PLOT LINES: Re(OC),Im(OC); Re(OA),Im(OA) !辺CA


LET a=ABS(BC) !辺BCの長さ
LET b=ABS(CA) !辺CA
LET c=ABS(AB) !辺AB

LET t=(a+b+c)/2 !ヘロンの公式より、面積S
LET S=SQR(t*(t-a)*(t-b)*(t-c))


!------------------------------

LET OI=(a*OA+b*OB+c*OC)/(a+b+c) !内心の位置ベクトル
SET AREA COLOR 4
DRAW disk WITH SCALE(0.2)*SHIFT(OI)

SET LINE COLOR 4
LET R=S/t !S=△IAB+△IBC+△ICA=1/2*c*r+1/2*a*r+1/2*b*r=r*(a+b+c)/2より
DRAW circle WITH SCALE(R)*SHIFT(OI) !内接円を描く　|z-α|=r

CALL DrawLine4(OA,OB,OC)
SUB DrawLine4(OA,OB,OC) !④∠BACを二等分する線
   FOR t=0 TO 10 !t=[0,∞]　※範囲は調整が必要である
      LET wAB=OB-OA
      LET wAC=OC-OA
      LET OP=OA+(fnNormalize(wAB)+fnNormalize(wAC))*t !ひし形AbDcの対角線AD
      PLOT LINES: Re(OP),Im(OP);
   NEXT t
   PLOT LINES
END SUB
CALL DrawLine4(OB,OC,OA)
CALL DrawLine4(OC,OA,OB)


!------------------------------

LET OG=(OA+OB+OC)/3 !重心の位置ベクトル
SET AREA COLOR 2
DRAW disk WITH SCALE(0.2)*SHIFT(OG)

SET LINE COLOR 2
CALL DrawLine2(OA,OB,OC)
SUB DrawLine2(OA,OB,OC) !②点Aと辺BCの中点を通る線
   FOR t=0 TO 10 !t=[0,∞]　※範囲は調整が必要である
      LET OP=OA+(OB+OC-2*OA)*t !平行四辺形ABDCの対角線AD
      PLOT LINES: Re(OP),Im(OP);
   NEXT t
   PLOT LINES
END SUB
CALL DrawLine2(OB,OC,OA)
CALL DrawLine2(OC,OA,OB)



!------------------------------

LET sin2A=2*S*(b^2+c^2-a^2)/(b*c)^2 !sin2A=2*cosA*sinA、余弦定理cosA=(b^2+c^2-a^2)/(2*b*c)、面積S=1/2*b*c*sinAより
LET sin2B=2*S*(c^2+a^2-b^2)/(c*a)^2
LET sin2C=2*S*(a^2+b^2-c^2)/(a*b)^2
LET OQ=(sin2A*OA+sin2B*OB+sin2C*OC)/(sin2A+sin2B+sin2C) !外心の位置ベクトル
SET AREA COLOR 3
DRAW disk WITH SCALE(0.2)*SHIFT(OQ)

SET LINE COLOR 3
LET R=a*b*c/(4*S) !正弦定理a/sinA=b/sinB=c/sinC=2*Rと面積S=1/2*b*c*sinAより
DRAW circle WITH SCALE(R)*SHIFT(OQ) !外接円を描く


!------------------------------

LET tanA=2*S/(b^2+c^2-a^2) !tanA=sinA/cosAと余弦定理cosA=(b^2+c^2-a^2)/(2*b*c)と面積S=1/2*b*c*sinAより
LET tanB=2*S/(c^2+a^2-b^2)
LET tanC=2*S/(a^2+b^2-c^2)
LET OH=(tanA*OA+tanB*OB+tanC*OC)/(tanA+tanB+tanC) !垂心の位置ベクトル
SET AREA COLOR 1
DRAW disk WITH SCALE(0.2)*SHIFT(OH)


FOR j=0 TO h !画面全体を走査する
   LET y=worldy(j) !ドットをxy座標に変換する
   FOR i=0 TO w
      LET x=worldx(i)

      LET OP=v(x,y)

      SET POINT COLOR 3 !外心
      IF ABS( fnDOT(OP-(OB+OC)/2,BC) )<cEps THEN PLOT POINTS: x,y !③辺BCの垂直二等分線
      IF ABS( fnDOT(OP-(OC+OA)/2,CA) )<cEps THEN PLOT POINTS: x,y
      IF ABS( fnDOT(OP-(OA+OB)/2,AB) )<cEps THEN PLOT POINTS: x,y

      SET POINT COLOR 1 !垂心
      IF ABS( fnDOT(OP-OA,BC) )<cEps THEN PLOT POINTS: x,y !①点Aから辺BCへの垂線　AP⊥BC
      IF ABS( fnDOT(OP-OB,CA) )<cEps THEN PLOT POINTS: x,y
      IF ABS( fnDOT(OP-OC,AB) )<cEps THEN PLOT POINTS: x,y

   NEXT i
NEXT j


END

LET t0=TIME


PUBLIC NUMERIC N !カードの枚数
LET N=13

DIM c(N) !最初のパターン　※
DATA 1,2,3,4,5,6,7,8,9,10,11,12,13
MAT READ c

PUBLIC NUMERIC GOAL(100) !最終のパターン　※
MAT GOAL=ZER(N)
DATA -9,-8,-7,1,2,4,5,6,-3,10,11,12,13
MAT READ GOAL

PUBLIC NUMERIC LIMIT !手数の上限　※
LET LIMIT=10

DIM A(LIMIT) !枚数
CALL backtrack(c,1,A)


PRINT "計算時間=";TIME-t0

END


EXTERNAL SUB reverse(c(),P) !先頭からP枚を裏返す
FOR i=1 TO INT(P/2) !交換
   LET t=c(i)
   LET c(i)=c(P-i+1)
   LET c(P-i+1)=t
NEXT i
FOR i=1 TO P !反転
   LET c(i)=-c(i)
NEXT i
END SUB

EXTERNAL SUB backtrack(c(),K,A()) !バックトラック法で検証する
IF K<=LIMIT THEN !手数の上限内なら
   DIM x(N)
   MAT x=c !save it

   FOR i=1 TO N !枚数を変える
      LET A(K)=i
      CALL reverse(c,i) !反転

      FOR j=1 TO N !最終のパターンかどうか確認する
         IF c(j)<>GOAL(j) THEN EXIT FOR
      NEXT j
      IF j>N THEN !一致したら
         LET LIMIT=K-1 !上限を狭める　※最初に見つかったもの
         PRINT K;"手"
         FOR j=1 TO K
            PRINT A(j); !枚数
         NEXT j
         PRINT
         !!!MAT PRINT c; !debug
      ELSE
         CALL backtrack(c,K+1,A) !次へ
      END IF

      MAT c=x !restore it
   NEXT i
END IF
END SUB

LET t0=TIME


PUBLIC NUMERIC N !カードの枚数
LET N=13

DIM c(N) !最初のパターン　※
DATA 1,2,3,4,5,6,7,8,9,10,11,12,13
MAT READ c

PUBLIC NUMERIC LIMIT !手数の上限　※
LET LIMIT=4

PUBLIC NUMERIC cMIN !最少手数
LET cMIN=LIMIT+1

DIM A(LIMIT) !枚数
CALL backtrack(c,0,A)

IF cMIN<=LIMIT THEN PRINT "最少手数=";cMIN


PRINT "計算時間=";TIME-t0

END


EXTERNAL SUB reverse(c(),P) !先頭からP枚を裏返す
FOR i=1 TO INT(P/2) !交換
   LET t=c(i)
   LET c(i)=c(P-i+1)
   LET c(P-i+1)=t
NEXT i
FOR i=1 TO P !反転
   LET c(i)=-c(i)
NEXT i
END SUB

EXTERNAL SUB backtrack(c(),K,A()) !バックトラック法で検証する
LET s=0
FOR j=1 TO N !表の枚数を確認する
   IF c(j)<0 THEN LET s=s+1
NEXT j
IF s=4 THEN !一致したら
   IF K<cMIN THEN LET cMIN=K !最少手数を記録する
   PRINT K;"手"
   FOR j=1 TO K
      PRINT A(j); !枚数列
   NEXT j
   PRINT
   MAT PRINT c; !最終のパターン
ELSE
   IF K<LIMIT THEN !手数の上限内なら
      DIM x(N)
      MAT x=c !save it

      FOR i=1 TO N !枚数を変える
         IF K>=1 AND i=A(K) THEN !同じ手が続く場合は無効！
         ELSE
            LET A(K+1)=i
            CALL reverse(c,i) !反転

            CALL backtrack(c,K+1,A) !次へ

            MAT c=x !restore it
         END IF
      NEXT i
   END IF
END IF
END SUB

!長方形を正方形で分割する
!例
!　┌──┬┬┐
!　│　　├┴┤
!　│　　│　│
!　└──┴─┘

LET a=65
LET b=47

IF a< b THEN !a≧bとする
   LET t=a !swap it
   LET a=b
   LET b=t
END IF

DO UNTIL b=0 !ユークリッドの互除法、連分数展開
   LET q=INT(a/b) !商
   PRINT "一辺";b;"の正方形が";q;"個"
   LET r=a-q*b !余り　※長方形をすべて正方形に分割できれば０となる
   LET a=b
   LET b=r
LOOP

PRINT "最大公約数=";a

END

一辺 47 の正方形が 1 個
一辺 18 の正方形が 2 個
一辺 11 の正方形が 1 個
一辺 7 の正方形が 1 個
一辺 4 の正方形が 1 個
一辺 3 の正方形が 1 個
一辺 1 の正方形が 3 個
最大公約数= 1

一辺 21 の正方形が 1 個
一辺 12 の正方形が 1 個
一辺 9 の正方形が 1 個
一辺 3 の正方形が 3 個

!長方形（正方形）の正方形分割

LET t0=TIME

PUBLIC NUMERIC a,b !長方形の縦、横の長さ
LET a=33
LET b=32

PUBLIC NUMERIC M !相異なる大きさの正方形
LET M=MIN(a,b)
DIM S(M) !未使用の個数
MAT S=CON


PRINT a;"×";b

DIM R(a,b) !元の長方形
MAT R=ZER

PUBLIC NUMERIC ANSWER_COUNT !解答数
LET ANSWER_COUNT=0

CALL backtrack(R,S,0,1,1)

IF ANSWER_COUNT=0 THEN PRINT "解なし"


PRINT "計算時間=";TIME-t0

END


EXTERNAL SUB backtrack(R(,),S(),W,x0,y0) !バックトラック法で検証する
!IF W>9 THEN EXIT SUB !個数の制限　※←←←←← ①

CALL serach(R,x0,y0,x,y) !配置位置
IF x>b AND y>a THEN !すべて埋め尽くされたら
   LET ANSWER_COUNT=ANSWER_COUNT+1
   PRINT "No.";ANSWER_COUNT; "個数=";W
   MAT PRINT USING(REPEAT$(" ##",b)): R

ELSE
   FOR i=M TO 1 STEP -1 !一辺ｉの正方形
      IF S(i)>0 THEN !未使用なら　※←←←←← ②

         CALL check(R,x,y,i, ok)
         IF ok=1 THEN !配置可能なら

            LET S(i)=S(i)-1 !使用中
            FOR dy=0 TO i-1 !置いてみる
               LET t=y+dy
               FOR dx=0 TO i-1
                  LET R(t,x+dx)=i
               NEXT dx
            NEXT dy

            CALL backtrack(R,S,W+1,x,y) !次へ

            LET S(i)=S(i)+1 !未使用
            FOR dy=0 TO i-1 !取り除く
               LET t=y+dy
               FOR dx=0 TO i-1
                  LET R(t,x+dx)=0
               NEXT dx
            NEXT dy

         END IF

      END IF !　※←←←←←
   NEXT i
END IF
END SUB

EXTERNAL SUB serach(R(,),x0,y0, x,y) !空き位置を探す
LET y=y0 !その行
FOR x=x0 TO b
   IF R(y0,x)=0 THEN EXIT SUB
NEXT x
FOR y=y0+1 TO a !以降の行
   FOR x=1 TO b
      IF R(y,x)=0 THEN EXIT SUB
   NEXT x
NEXT y
END SUB

EXTERNAL SUB check(R(,),x,y,c, ok) !位置(x,y)に一辺cの正方形が配置可能かどうか確認する
LET ok=0
IF x+c-1>b THEN EXIT SUB !はみ出す
IF y+c-1>a THEN EXIT SUB
FOR dy=0 TO c-1
   LET t=y+dy
   FOR dx=0 TO c-1
      IF R(t,x+dx)<>0 THEN EXIT SUB !十分な空きがない
   NEXT dx
NEXT dy
LET ok=1 !可能！
END SUB

!問題
!　ゲーム「テトリス」のテトロミノの回転の実装
!　例　Ｔ型
!　　　■
!　　■■　→　■■■
!　　　■　　　　■

!答え
!PxQ行列の成分を定義する

LET P=3 !行
LET Q=2 !列

DIM S$(P*Q) !パターン
DATA "　","■"
DATA "■","■"
DATA "　","■"
MAT READ S$

!PxQ行列の成分を90°回転させる場合、成分の番号付けは
!　 1 2　　　 2 4 6
!　 3 4　→　 1 3 5
!　 5 6
!となる。

!成分の並びを入れ替える置換σ、写像fを定義する

DIM T(Q,P) !行列の成分番号
FOR i=1 TO Q !反時計まわりに90°回転
   FOR j=1 TO P
      LET T(i,j)=Q*(j-1)+(Q-i+1)
   NEXT j
NEXT i

CALL PrintOutA(S$,T,Q,P)
SUB PrintOutA(T$(),A(,),P,Q) !並び順AによるP行Q列の行列を表示する
!MAT PRINT A; !成分の番号
   FOR i=1 TO P
      FOR j=1 TO Q
         PRINT S$(A(i,j));
      NEXT j
      PRINT
   NEXT i
END SUB



!次に個々の並べ替えをみていくことにする。

LET P=5 !行
LET Q=3 !列

DIM B(P,Q),tB(Q,P) !行列の成分番号
FOR i=1 TO P
   FOR j=1 TO Q
      LET B(i,j)=Q*(i-1)+j !元の成分(i,j)
   NEXT j
NEXT i
MAT PRINT B;

FOR i=1 TO Q !転置
   FOR j=1 TO P
      LET tB(i,j)=Q*MOD(j-1,P)+i !(MOD(j-1,P)+1,i)
   NEXT j
NEXT i
MAT PRINT tB;


FOR i=1 TO Q !反時計まわりに90°回転
   FOR j=1 TO P
      LET tB(i,j)=Q*(j-1)+(Q-i+1) !(j,Q-i+1)
   NEXT j
NEXT i
MAT PRINT tB;

FOR i=1 TO P !180°
   FOR j=1 TO Q
      LET B(i,j)=Q*(P-i)+(Q-j+1) !(P-i+1,Q-j+1)
   NEXT j
NEXT i
MAT PRINT B;

FOR i=1 TO Q !270°,-90°
   FOR j=1 TO P
      LET tB(i,j)=Q*(P-j)+i !(P-j+1,i)
   NEXT j
NEXT i
MAT PRINT tB;


FOR i=1 TO P !Ｘ軸対称
   FOR j=1 TO Q
      LET B(i,j)=Q*(P-i)+j !(P-i+1,j)
   NEXT j
NEXT i
MAT PRINT B;

FOR i=1 TO P !Ｙ軸対称
   FOR j=1 TO Q
      LET B(i,j)=Q*(i-1)+(Q-j+1) !(i,Q-j+1)
   NEXT j
NEXT i
MAT PRINT B;


END

!「同じものを含む順列」に番号をつける方法
!例 aabc,aacb,abac,…,caba,cbaa　(2+1+1)!/(2!*1!*1!)=12通り


LET M=3 !異なるｍ種類
DIM B(M) !それぞれの個数
DATA 2,2,1 !※{1,1,2,2,3}の意
MAT READ B

LET N=0 !総数
FOR i=1 TO M
   LET N=N+B(i)
NEXT i


DIM A(N)
DATA 1,1,2,2,3 !最初のパターン
!!!DATA 3,2,2,1,1 !最後のパターン　※前の順列　←←←←←
MAT READ A

FOR i=0 TO PermFactorialM(B,M)-1
!!!FOR i=PermFactorialM(B,M)-1 TO 0 STEP -1 !　※前の順列　←←←←←

   PRINT "No.";i
   MAT PRINT A;
   PRINT PermFactorialM2Num(A,N,B,M) !符号化

   DIM AA(N)
   MAT AA=ZER
   CALL Num2PermFactorialM(i, AA,N,B,M) !復号化
   MAT PRINT AA;

   CALL NextPermFactorial(A,N, rc) !次へ
   PRINT "rc=";rc

NEXT i


END


!●同じものをそれぞれp,q,…,r個ずつ含む総数n個（n=p+q+ … +r）のものをすべて並べる

!(p+q+ … +r)!/(p!*q!* … *r!)通りの順列パターン ⇔ 0～(p+q+ … +r)!/(p!*q!* … *r!)の番号

EXTERNAL FUNCTION PermFactorialM(B(),M) !同じものを含む順列の「場合の数」
LET s=B(1) !総数 p+q+ … +r
LET t=FACT(B(1)) !階乗の積 p!*q!* … *r!
FOR i=2 TO M
   LET s=s+B(i)
   LET t=t*FACT(B(i))
NEXT i
LET PermFactorialM=FACT(s)/t
END FUNCTION


EXTERNAL FUNCTION PermFactorialM2Num(A(),N,B(),M) !順列パターンに番号を付ける　※辞書式順序
DIM w(M)
MAT w=B !copy it
LET v=0
FOR i=1 TO N-1 !左端に着目する
   LET t=A(i)
   FOR j=1 TO t-1 !左端を１～A(i)-1として
      IF w(j)>0 THEN !その左端を除いた残りで最大の順列をつくる
         LET w(j)=w(j)-1
         !!!MAT PRINT w; !debug
         LET v=v+PermFactorialM(w,M) !その順列の番号を求める
         LET w(j)=w(j)+1
      END IF
   NEXT j
   LET w(t)=w(t)-1 !左端を消して、次へ
NEXT i
LET PermFactorialM2Num=v
END FUNCTION

EXTERNAL SUB Num2PermFactorialM(h, A(),N,B(),M) !番号から順列パターンを生成する　※辞書式順序
DIM w(M)
MAT w=B !copy it
LET v=h
FOR i=1 TO N
   FOR j=1 TO M !左端を１～Mとして
      IF w(j)>0 THEN !その左端を除いた残りで最大の順列をつくる
         LET w(j)=w(j)-1
         !!!MAT PRINT w;
         LET t=PermFactorialM(w,M) !その順列の番号を求める
         IF v<t THEN EXIT FOR
         LET v=v-t
         LET w(j)=w(j)+1
      END IF
   NEXT j
   LET A(i)=j !次へ
NEXT i
END SUB


EXTERNAL SUB NextPermFactorial(A(),N, rc) !辞書式順序で次の順列を返す　※「異なるｎ個のもの」と共通
LET i=N-1 !順列を右から左にみて、増加列から減少列に変わる位置iを探す
DO WHILE i>0 AND A(i)>=A(i+1) !０は番人
!!!DO WHILE i>0 AND A(i)<=A(i+1) !０は番人　※前の順列　←←←←←
   LET i=i-1
LOOP
IF i=0 THEN !A(1)>A(2)>A(3)> … >A(N)なら　例. N=4、4,3,2,1
   LET rc=0 !完了
   EXIT SUB
END IF

LET j=N !その位置iより右で、A(i)以上で最小の数A(j)を探す
DO WHILE A(i)>=A(j)
!!!DO WHILE A(i)<=A(j) !　※前の順列　←←←←←
   LET j=j-1
LOOP
LET t=A(i) !A(i)とA(j)を交換する
LET A(i)=A(j)
LET A(j)=t

LET i=i+1 !(i+1)からNまでの範囲を逆順にする
LET j=N
DO WHILE i<j
   LET t=A(i) !swap it
   LET A(i)=A(j)
   LET A(j)=t
   LET i=i+1
   LET j=j-1
LOOP
LET rc=1 !未了
END SUB

!制約式
!　a11*x1+a12*x2+ … +a1n*xn≧b1　※≦
!　a21*x1+a22*x2+ … +a2n*xn≧b2　※≦
!　　　：
!　am1*x1+am2*x2+ … +amn*xn≧bm　※≦
!目的関数
!　z=c1*x1+c2*x2+ … *cn*xn

!例
!　　x1+3*x2≦18
!　2*x1+3*x2≦21
!　3*x1+　x2≦21
!　　x1,x2≧0
!　z=2*x1+x2の最大値

LET N=2 !変数の数
LET M=3 !制約式の数

LET R=M+1 !行数
LET C=M+N+1 !桁数

DIM A(R,C) !係数行列
DATA  1, 3, 1,0,0, 18
DATA  2, 3, 0,1,0, 21
DATA  3, 1, 0,0,1, 21
DATA -2,-1, 0,0,0,  0
MAT READ A

DO
   LET mn=9999 !列選択
   FOR k=1 TO C-1
      IF A(R,k)<mn THEN
         LET mn=A(R,k)
         LET y=k
      END IF
   NEXT k
   IF mn>=0 THEN EXIT DO
   LET mn=9999 !行選択
   FOR k=1 TO R-1
      LET p=A(k,C)/A(k,y)
      IF A(k,y)>0 AND p<mn THEN
         LET mn=p
         LET x=k
      END IF
   NEXT k
   LET p=A(x,y) !ピボット係数
   FOR k=1 TO C !ピボット行をｐで割る
      LET A(x,k)=A(x,k)/p
   NEXT k
   FOR k=1 TO R !ピボット列の掃き出し
      IF k<>x THEN
         LET d=A(k,y)
         FOR j=1 TO C
            LET A(k,j)=A(k,j)-d*A(x,j)
         NEXT j
      END IF
   NEXT k

LOOP
MAT PRINT A; !debug

FOR k=1 TO N
   LET flg=-1
   FOR j=1 TO R
      IF A(j,k)=1 THEN LET flg=j
   NEXT j
   IF flg<>-1 THEN PRINT "x";STR$(k);"=";A(flg,C)
NEXT k
PRINT "z=";A(R,C)

END

DEF f(x1,x2)=x1^2-x1*x2+3*x2^2 !最小化する関数

DEF df1(x1,x2)=2*x1-x2 !∂f/∂x1
DEF df2(x1,x2)=-x1+6*x2 !∂f/∂x2

LET N=2 !変数の数

DIM H(N,N) !ヘッセ行列　※正則かつ正定値行列（固有値がすべて正の値）
DATA  2,-1
DATA -1, 6
MAT READ H

DIM x(N) !初期値
DATA 1,2
MAT READ x

LET cEPS=1E-8 !精度

DIM iH(N,N)
MAT iH=INV(H)

LET iter=200 !繰り返し回数
FOR i=1 TO iter
   PRINT "x1=";x(1), "x2=";x(2), "f=";f(x(1),x(2))

   DIM df(N) !勾配∇f
   LET df(1)=df1(x(1),x(2))
   LET df(2)=df2(x(1),x(2))
   !!!MAT PRINT df; !debug

   DIM T(N)
   MAT T=iH*df !変位Δx
   MAT T=x-T !近似解x[k+1]=x[k]-Δx
   !!!MAT PRINT T; !debug

   DIM w(N) !近似解x[k+1]とx[k]のノルム
   MAT w=T-x
   IF DOT(w,w)<cEPS*cEPS THEN EXIT FOR !収束したなら、終了！

   MAT x=T !次へ
NEXT i

END

!１～ｎの数が記入されたカードがｍ枚ずつある。
!これを１～ｍ*ｎの位置に置くとき、その配置の「場合の数」を求めよ。
!ただし、番号ｎのカードどうしはｎだけ離れて置く必要がある。
!例　ｎ=４、ｍ=３で、４のカードの場合
!　123456789ABC←位置
!　4++++4++++4+
!　+4++++4++++4　２通り

!答え
!m=3
!n=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15, …
!　0,0,0,0,0,0,0,0,6,10, 0, 0, 0,??,??, …
!※パソコンの性能により調査できず

!PentiumⅢ700MHz,WindowsMe,２進モード　n=13の場合、計算時間= 175.71 秒

LET t0=TIME

PUBLIC NUMERIC N,M
LET N=9 !１～ｎのカード
LET M=3 !ｍ枚ずつ

PUBLIC NUMERIC ANSWER_COUNT
LET ANSWER_COUNT=0

DIM A(M*N) !配置位置
MAT A=ZER
CALL backtrack(N,A) !N,…,2,1　計算量O(FACT(N))

IF ANSWER_COUNT=0 THEN PRINT "解なし"

PRINT "計算時間=";TIME-t0

END

EXTERNAL SUB backtrack(p,A()) !バックトラック法で検証する
FOR i=1 TO M*N-(p+1)*(M-1) !カードpの１枚目を位置iに置く
   FOR k=0 TO M-1
      IF A(i+(p+1)*k)<>0 THEN EXIT FOR
   NEXT k
   IF k=M THEN !ｍ枚とも置けるなら

      FOR k=0 TO M-1 !仮に置いてみる
         LET A(i+(p+1)*k)=p
      NEXT k

      IF p=1 THEN !すべて置けたら、結果を表示する
         LET ANSWER_COUNT=ANSWER_COUNT+1 !解答数
         PRINT "No.";ANSWER_COUNT
         MAT PRINT A; !debug
      ELSE
         CALL backtrack(p-1,A) !次へ
      END IF

      FOR k=0 TO M-1 !元に戻す
         LET A(i+(p+1)*k)=0
      NEXT k

   END IF
NEXT i
END SUB

No. 1
 1  9  1  6  1  8  2  5  7  2  6  9  2  5  8  4  7  6  3  5  4  9  3  8  7  4  3

No. 2
 1  9  1  2  1  8  2  4  6  2  7  9  4  5  8  6  3  4  7  5  3  9  6  8  3  5  7

No. 3
 1  8  1  9  1  5  2  6  7  2  8  5  2  9  6  4  7  5  3  8  4  6  3  9  7  4  3

No. 4
 3  4  7  9  3  6  4  8  3  5  7  4  6  9  2  5  8  2  7  6  2  5  1  9  1  8  1

No. 5
 7  5  3  8  6  9  3  5  7  4  3  6  8  5  4  9  7  2  6  4  2  8  1  2  1  9  1

No. 6
 3  4  7  8  3  9  4  5  3  6  7  4  8  5  2  9  6  2  7  5  2  8  1  6  1  9  1

DIM A(10)
LET A(5)=5
MAT A=(A(5))*CON
MAT PRINT A;

DIM B(10,10)
LET B(5,5)=5
MAT B=(B(5,5))*IDN
MAT PRINT B;

END

!非線形関数 f(x1,x2)=(x1-2)^2+(x1-2*x2)^2 の最小値
!初期値を(x1,x2)=(0,0)、ステップ幅α=0.1 とする。

DEF f(x1,x2)=(x1-2)^2+(x1-2*x2)^2 !関数

DEF df1(x1,x2)=2*(x1-2)+2*(x1-2*x2) !∂f/∂x1
DEF df2(x1,x2)=-4*(x1-2*x2) !∂f/∂x2

LET N=2 !変数の数

DIM x(N) !初期値
DATA 1,2
MAT READ x

LET a=0.1 !ステップ幅

LET cEPS=1E-8 !精度

LET iter=200 !繰り返し回数
FOR k=1 TO iter
   PRINT "x1=";x(1), "x2=";x(2), "f=";f(x(1),x(2))

   DIM df(N) !勾配∇f
   LET df(1)=df1(x(1),x(2))
   LET df(2)=df2(x(1),x(2))
   !!!MAT PRINT df; !debug

   LET norm=DOT(df,df)
   IF norm<cEPS*cEPS THEN EXIT FOR !勾配が０なら、終了！

   DIM S(N)
   MAT S=(-1)*df !降下方向 s[k]=-∇f(x[k])

   DIM w(N),s0(N)
   IF k<>1 THEN
   !λ[k]=∥∇f(x[k])∥^2/∥∇f(x[k-1]∥^2　フレッチャー・リーブス（Fletcher-Reeves）の公式
      LET l=norm/norm0

      MAT w=l*s0 !s[k]=-∇f(x[k])+λ[k-1]*s[k-1]
      MAT S=S+w
   END IF
   LET norm0=norm
   MAT s0=S

   MAT w=a*S !近似解 x[k+1]=x[k]+α*s[k]
   MAT x=x+w
NEXT k

END

!PERT図（アローダイアグラム）

SET WINDOW -0.5,6,-0.5,6 !作画領域
DRAW grid
SET TEXT JUSTIFY "center","half"

LET N=9 !結合点（ノード）の数

DATA 1,3 !結合点1のx,y座標
DATA 2,4 !2
DATA 2,3 !3
DATA 2,1 !4
DATA 3,1 !5
DATA 3,3 !6
DATA 3,4 !7
DATA 4,3 !8
DATA 5,3 !9

DIM xNd(N),yNd(N) !結合点のx,y座標
FOR i=1 TO N
   READ xNd(i),yNd(i)
NEXT i


LET T=11 !工程の数
DATA 1,2, 5,"A" !工程Aの結合点と日数　※始端:1、終端:N
DATA 2,7, 2,"B"
DATA 7,8, 1,"D"
DATA 1,3,15,"E"
DATA 3,6, 7,"F"
DATA 6,8, 7,"H"
DATA 1,4,20,"C"
DATA 4,5, 1,"G"
DATA 5,6, 0,""
DATA 5,8, 3,"i"
DATA 8,9, 1,"J"

DIM M(N,N) !有向グラフの隣接行列　※成分(i,j)は、上三角行列:i→j、下三角行列:j→iの意
MAT M=(-1)*CON
FOR i=1 TO T
   READ n1,n2,d,nm$
   IF n1<1 OR n2<1 OR n1>N OR n2>N THEN
      PRINT "結合点番号が１～";N;"の範囲ではありません。", n1;n2;d;nm$
      STOP
   END IF
   IF d<0 THEN
      PRINT "日数が負です。", n1;n2;d;nm$
      STOP
   END IF
   LET M(n1,n2)=d
   LET M(n2,n1)=d

   IF d<>0 THEN SET LINE STYLE 1 ELSE SET LINE STYLE 3 !実線、破線
   PLOT LINES: xNd(n1),yNd(n1); xNd(n2),yNd(n2) !結線
   PLOT TEXT ,AT (xNd(n1)+xNd(n2))/2,(yNd(n1)+yNd(n2))/2+0.15: nm$&"("&STR$(d)&")" !工程（日数）
NEXT i
MAT PRINT M;


PUBLIC NUMERIC Z(100),ZZ(100) !結合点（ノード）の日数

!最早結合点時刻
DIM R(N),A(N) !経路、日数
MAT A=ZER
MAT R=ZER
LET R(1)=1 !始端１から開始する
MAT Z=ZER(N)
CALL backtrack(1,M,N,A,R,0)
MAT PRINT Z; !結果

PRINT


!最遅結合点時刻
LET mx=Z(N) !終端での日数
MAT A=(mx)*CON
MAT R=ZER
LET R(N)=N !終端Ｎから開始する
MAT ZZ=A
CALL backtrack2(N,M,N,A,R,(A(N))) !※A(N)は値渡し
MAT PRINT ZZ; !結果


FOR i=1 TO N !結合点を描く
   PLOT TEXT ,AT xNd(i),yNd(i)-0.25: STR$(Z(i))&"/"&STR$(ZZ(i)) !日数

   IF Z(i)=ZZ(i) THEN SET AREA COLOR 4 ELSE SET AREA COLOR 1 !クリティカルパス
   DRAW disk WITH SCALE(0.1)*SHIFT(xNd(i),yNd(i))

   PLOT TEXT, AT xNd(i),yNd(i)+0.25: STR$(i) !番号
NEXT i

END


EXTERNAL SUB backtrack(c,M(,),N,A(),R(),s) !バックトラック法で全経路を調査する
IF c=N THEN !終端Ｎなら
   MAT PRINT R; !debug
   MAT PRINT A;
   FOR i=1 TO N !各ノードの日数を記録する
      IF A(i)>Z(i) THEN LET Z(i)=A(i) !最大のもの
   NEXT i
ELSE
   FOR j=c TO N !次の経路を探す　※上三角行列
      LET w=M(c,j)
      IF w>=0 THEN !経路c→jがあれば
         LET s=s+w !累計
         LET R(j)=j !経路を記録する
         LET A(j)=s !日数
         CALL backtrack(j,M,N,A,R,s) !次へ
         LET s=s-w !元に戻す
         LET R(j)=0
         LET A(j)=A(1)
      END IF
   NEXT j
END IF
END SUB

EXTERNAL SUB backtrack2(c,M(,),N,A(),R(),s) !バックトラック法で全経路を調査する
IF c=1 THEN !始端１なら
   MAT PRINT R; !debug
   MAT PRINT A;
   FOR i=1 TO N !各ノードの日数を記録する
      IF A(i)<ZZ(i) THEN LET ZZ(i)=A(i) !最小のもの
   NEXT i
ELSE
   FOR j=1 TO c !次の経路を探す　※下三角行列
      LET w=M(c,j)
      IF w>=0 THEN !経路j→cがあれば
         LET s=s-w !累計
         LET R(j)=j !経路を記録する
         LET A(j)=s !日数
         CALL backtrack2(j,M,N,A,R,s) !次へ
         LET s=s+w !元に戻す
         LET R(j)=0
         LET A(j)=A(N)
      END IF
   NEXT j
END IF
END SUB

OPEN #1: NAME "TEST.TXT" ,ACCESS OUTIN
ERASE #1
LET a=7
LET b=15
PRINT #1: a, ","; b !ケース１
!!!PRINT #1: a; ","; b !ケース２
CLOSE #1
END

 7                      , 15

 7 , 15

OPEN #1: NAME "TEST.TXT" ,ACCESS INPUT
INPUT #1: a,b
PRINT a,b
CLOSE #1
END

順列

●異なるn個のものをすべて並べる　※{1,2,3,…,n}上の全順列（n-順列）

FACT(N)　場合の数　※組み込み関数
FUNCTION PermFactorial2Num(A(),N)　パターンに番号を付ける
SUB Num2PermFactorial(h, A(),N)　番号からパターンを生成する
SUB PrevPermFactorial(A(),N, rc)　辞書式順序で１つ前を返す
SUB NextPermFactorial(A(),N, rc)　辞書式順序で１つ次を返す


●異なるn個のものから重複を許さずr個を取り出して並べる　※{1,2,3,…,n}上のr-順列

PERM(N,R)　※組み込み関数
FUNCTION Perm2Num(A(),N,R)
SUB Num2Perm(h, A(),N,R)
SUB PrevPerm(A(),N,R, rc)
SUB NextPerm(A(),N,R, rc)


●同じものをそれぞれp,q,…,r個ずつ含む総数n個（n=p+q+ … +r）のものをすべて並べる

FUNCTION PermFactorialM(B(),M)　※(p+q+ … +r)!/(p!*q!* … *r!)
FUNCTION PermFactorialM2Num(A(),N,B(),M)
SUB Num2PermFactorialM(h, A(),N,B(),M)
SUB PrevPermFactorial(A(),N, rc)　※全順列と同じ
SUB NextPermFactorial(A(),N, rc)　※全順列と同じ


●異なるn個のものから重複を許してr個を取り出して並べる（重複順列）

N^R
FUNCTION ReptPerm2Num(A(),N,R)
SUB Num2ReptPerm(h, A(),N,R)
SUB PrevReptPerm(A(),N,R, rc)
SUB NextReptPerm(A(),N,R, rc)


●異なるn個のものをすべて円形に並べる（円順列）

FACT(N-1)　※組み込み関数
？
？
？
？


●異なるn個のものをすべて円形に並べる（左右対称なものは除く）（数珠順列）

FACT(N-1)/2　※組み込み関数
？
？
？
？


組合せ

●異なるn個のものから重複を許さずr個を取り出す組合せ

COMB(N,R)　※組み込み関数
FUNCTION Comb2Num(A(),N,R)
SUB Num2Comb(h, A(),N,R)
SUB PrevComb(A(),N,R, rc)
SUB NextComb(A(),N,R, rc)


●異なるn個のものから重複を許してr個を取り出す組合せ（重複組合せ）

COMB(N+R-1,R)　※組み込み関数
FUNCTION ReptComb2Num(A(),N,R)
SUB Num2ReptComb(h, A(),N,R)
SUB PrevReptComb(A(),N,R, rc)
SUB NextReptComb(A(),N,R, rc)

LET N=5
LET R=3

DIM A(R) !最初のパターン
FOR i=1 TO R !初期値 A={1,1,1,…,1}　※重複組合せ
   LET A(i)=1
NEXT i

LET s=0 !番号付け　０～
DO
   PRINT "No.";s
   MAT PRINT A;

   PRINT ReptComb2Num(A,N,R) !符号化

   DIM AA(R)
   MAT AA=ZER
   CALL Num2ReptComb(s,AA,N,R) !復号化
   MAT PRINT AA;

   CALL NextReptComb(A,N,R, rc) !次へ
   LET s=s+1
LOOP WHILE rc<>0
PRINT


LET s=COMB(N+R-1,R)-1 !場合の数　nHr=n+r-1Cr
DO
   PRINT "No.";s
   MAT PRINT A;
   CALL PrevReptComb(A,N,R, rc)
   LET s=s-1
LOOP WHILE rc<>0


END


!最小完全ハッシュ関数 factoradic

!●異なるn個のものをすべて並べる　※{1,2,3,…,n}上の全順列（n-順列）

!n!通りの順列パターン ⇔ 0～(n!-1)の番号

EXTERNAL FUNCTION PermFactorial2Num(A(),N) !順列パターンに番号を付ける　※辞書式順序
FOR j=1 TO N-1 !階乗進数の各桁の値＋１　A[1..N]=(N-1)! … 3! 2! 1! 0!
   FOR k=j+1 TO N
      IF A(k)>=A(j) THEN LET A(k)=A(k)-1
   NEXT k
NEXT j
LET v=0
FOR j=N TO 1 STEP -1 !非負の10進数整数へ
   LET v=v*j+A(N-j+1)-1
NEXT j
LET PermFactorial2Num=v
END FUNCTION

EXTERNAL SUB Num2PermFactorial(h, A(),N) !番号から順列パターンを生成する　※辞書式順序
LET v=h !非負の10進数整数を階乗進数へ
FOR j=1 TO N
   LET t=INT(v/j)
   LET A(N-j+1)=v-t*j +1 !階乗進数の各桁の値＋１　A[1..N]=(N-1)! … 3! 2! 1! 0!
   LET v=t
NEXT j
FOR j=N-1 TO 1 STEP -1 !順列パターンへ
   FOR k=j+1 TO N
      IF A(k)>=A(j) THEN LET A(k)=A(k)+1
   NEXT k
NEXT j
END SUB


EXTERNAL SUB PrevPermFactorial(A(),N, rc) !辞書式順序で次の順列を返す
LET rc=0 !完了

LET i=N-1 !順列を右から左にみて、増加列から減少列に変わる位置iを探す
DO WHILE i>0 AND A(i)<=A(i+1) !０は番人
   LET i=i-1
LOOP
IF i=0 THEN EXIT SUB !A(1)>A(2)>A(3)> … >A(N)なら　例. N=4、4,3,2,1

LET j=N !その位置iより右で、A(i)以上で最小の数A(j)を探す
DO WHILE A(i)<=A(j)
   LET j=j-1
LOOP
LET t=A(i) !A(i)とA(j)を交換する
LET A(i)=A(j)
LET A(j)=t

LET i=i+1 !(i+1)からNまでの範囲を逆順にする
LET j=N
DO WHILE i<j
   LET t=A(i) !swap it
   LET A(i)=A(j)
   LET A(j)=t
   LET i=i+1
   LET j=j-1
LOOP
LET rc=1 !未了
END SUB

EXTERNAL SUB NextPermFactorial(A(),N, rc) !辞書式順序で次の順列を返す
LET rc=0 !完了

LET i=N-1 !順列を右から左にみて、増加列から減少列に変わる位置iを探す
DO WHILE i>0 AND A(i)>=A(i+1) !０は番人
!!!DO WHILE i>0 AND A(i)<=A(i+1) !０は番人　※前の順列　←←←←←
   LET i=i-1
LOOP
IF i=0 THEN EXIT SUB !A(1)>A(2)>A(3)> … >A(N)なら　例. N=4、4,3,2,1

LET j=N !その位置iより右で、A(i)以上で最小の数A(j)を探す
DO WHILE A(i)>=A(j)
!!!DO WHILE A(i)<=A(j) !　※前の順列　←←←←←
   LET j=j-1
LOOP
LET t=A(i) !A(i)とA(j)を交換する
LET A(i)=A(j)
LET A(j)=t

LET i=i+1 !(i+1)からNまでの範囲を逆順にする
LET j=N
DO WHILE i<j
   LET t=A(i) !swap it
   LET A(i)=A(j)
   LET A(j)=t
   LET i=i+1
   LET j=j-1
LOOP
LET rc=1 !未了
END SUB


!●異なるn個のものから重複を許さずr個を取り出して並べる　※{1,2,3,…,n}上のr-順列

!perm(n,r)通りの順列パターン ⇔ 0 ～ perm(n,r)-1 の番号

EXTERNAL FUNCTION Perm2Num(A(),N,R) !順列パターンに番号を付ける　※辞書式順序
FOR j=1 TO R-1 !(階乗)進数の各桁の値＋１　A[1..R]=PERM(N-1,R-1) … PERM(N-j,R-j) … PERM(N-R,0)
   FOR k=j+1 TO R
      IF A(k)>=A(j) THEN LET A(k)=A(k)-1
   NEXT k
NEXT j
LET v=0
FOR j=R TO 1 STEP -1 !非負の10進数整数へ
   LET v=v*(N-R+j)+A(R-j+1)-1
NEXT j
LET Perm2Num=v
END FUNCTION

EXTERNAL SUB Num2Perm(h, A(),N,R) !番号から順列パターンを生成する　※辞書式順序
LET v=h !非負の10進数整数を(階乗)進数へ
FOR j=1 TO R
   LET w=N-R+j
   LET t=INT(v/w)
   LET A(R-j+1)=v-t*w +1 !(階乗)進数の各桁の値＋１　A[1..R]=PERM(N-1,R-1) … PERM(N-j,R-j) … PERM(N-R,0)
   LET v=t
NEXT j
FOR j=R-1 TO 1 STEP -1 !順列パターンへ
   FOR k=j+1 TO R
      IF A(k)>=A(j) THEN LET A(k)=A(k)+1
   NEXT k
NEXT j
END SUB


EXTERNAL SUB PrevPerm(A(),N,R, rc) !辞書式順序で前の順列を返す
LET rc=0 !完了

FOR j=1 TO R-1 !(階乗)進数へ
   FOR k=j+1 TO R
      IF A(k)>=A(j) THEN LET A(k)=A(k)-1
   NEXT k
NEXT j

FOR i=R TO 1 STEP -1 !(階乗)進数で－１する
   IF A(i)>1 THEN !1～(N-i+1)で更新する
      LET A(i)=A(i)-1
      FOR k=i+1 TO R !ｉ桁より下の位を「最後の並び」にする
         LET A(k)=N-k+1
      NEXT k

      LET rc=1 !未了
      EXIT FOR
   END IF
NEXT i

FOR j=R-1 TO 1 STEP -1 !順列パターンへ
   FOR k=j+1 TO R
      IF A(k)>=A(j) THEN LET A(k)=A(k)+1
   NEXT k
NEXT j
END SUB

EXTERNAL SUB NextPerm(A(),N,R, rc) !辞書式順序で次の順列を返す
LET rc=0 !完了

FOR j=1 TO R-1 !(階乗)進数へ
   FOR k=j+1 TO R
      IF A(k)>=A(j) THEN LET A(k)=A(k)-1
   NEXT k
NEXT j

FOR i=R TO 1 STEP -1 !(階乗)進数で＋１する
   IF A(i)<N-i+1 THEN !1～(N-i+1)で更新する
      LET A(i)=A(i)+1
      FOR k=i+1 TO R !ｉ桁より下の位を「最初の並び」にする
         LET A(k)=1
      NEXT k

      LET rc=1 !未了
      EXIT FOR
   END IF
NEXT i

FOR j=R-1 TO 1 STEP -1 !順列パターンへ
   FOR k=j+1 TO R
      IF A(k)>=A(j) THEN LET A(k)=A(k)+1
   NEXT k
NEXT j
END SUB

!●同じものをそれぞれp,q,…,r個ずつ含む総数n個（n=p+q+ … +r）のものをすべて並べる

!(p+q+ … +r)!/(p!*q!* … *r!)通りの順列パターン ⇔ 0～(p+q+ … +r)!/(p!*q!* … *r!)-1の番号

EXTERNAL FUNCTION PermFactorialM(B(),M) !同じものを含む順列の「場合の数」
LET s=B(M) !総数 r, … ,q+ … +r,p+q+ … +r
LET t=1 !組合せ comb(r,r), … ,comb(q+ … +r,q),comb(p+q+ … +r,p)
FOR i=M-1 TO 1 STEP -1
   LET s=s+B(i)
   LET t=t*COMB(s,B(i)) !組合せ順列
NEXT i
LET PermFactorialM=t
!!!別解
!!!LET s=B(1) !総数 p+q+ … +r
!!!LET t=FACT(B(1)) !階乗の積 p!*q!* … *r!
!!!FOR i=2 TO M
!!!   LET s=s+B(i)
!!!   LET t=t*FACT(B(i))
!!!NEXT i
!!!LET PermFactorialM=FACT(s)/t
END FUNCTION


EXTERNAL FUNCTION PermFactorialM2Num(A(),N,B(),M) !順列パターンに番号を付ける　※辞書式順序
DIM w(M)
MAT w=B !copy it
LET v=0
FOR i=1 TO N-1 !左端に着目する
   LET t=A(i)
   FOR j=1 TO t-1 !左端を１～A(i)-1として
      IF w(j)>0 THEN !その左端を除いた残りで最大の順列をつくる
         LET w(j)=w(j)-1
         !!!MAT PRINT w; !debug
         LET v=v+PermFactorialM(w,M) !その順列の番号を求める
         LET w(j)=w(j)+1
      END IF
   NEXT j
   LET w(t)=w(t)-1 !左端を消して、次へ
NEXT i
LET PermFactorialM2Num=v
END FUNCTION

EXTERNAL SUB Num2PermFactorialM(h, A(),N,B(),M) !番号から順列パターンを生成する　※辞書式順序
DIM w(M)
MAT w=B !copy it
LET v=h
FOR i=1 TO N
   FOR j=1 TO M !左端を１～Mとして
      IF w(j)>0 THEN !その左端を除いた残りで最大の順列をつくる
         LET w(j)=w(j)-1
         !!!MAT PRINT w;
         LET t=PermFactorialM(w,M) !その順列の番号を求める
         IF v<t THEN EXIT FOR
         LET v=v-t
         LET w(j)=w(j)+1
      END IF
   NEXT j
   LET A(i)=j !次へ
NEXT i
END SUB


!●異なるn個のものから重複を許してr個を取り出して並べる（重複順列）

!N^R通りの順列パターン ⇔ 0 ～ N^R-1 の番号

EXTERNAL FUNCTION ReptPerm2Num(A(),N,R) !順列パターンに番号を付ける　※辞書式順序
LET v=A(1)-1
FOR i=2 TO R !N進法R桁の数を10進法へ
   LET v=v*N+A(i)-1
NEXT i
LET ReptPerm2Num=v
END FUNCTION

EXTERNAL SUB Num2ReptPerm(h, A(),N,R) !番号から順列パターンを生成する　※辞書式順序
LET v=h
LET i=R !桁位置
DO UNTIL v=0 !10進法の数をN進法R桁へ
   LET t=INT(v/N)
   LET A(i)=v-t*N+1 !1～N　※剰余
   LET v=t
   LET i=i-1
LOOP
FOR k=i TO 1 STEP -1 !残りの桁を埋める
   LET A(k)=1
NEXT k
END SUB


EXTERNAL SUB PrevReptPerm(A(),N,R, rc) !辞書式順序で前の順列を返す
LET rc=0 !完了
FOR i=R TO 1 STEP -1 !N進法R桁の数として
   IF A(i)>1 THEN !1～Nで更新する
      LET A(i)=A(i)-1
      FOR j=i+1 TO R !A(i)=A(i+1)= … =A(R)　最後の並び
         LET A(j)=N
      NEXT j
      LET rc=1 !未了
      EXIT SUB
   END IF
NEXT i
END SUB

EXTERNAL SUB NextReptPerm(A(),N,R, rc) !辞書式順序で次の順列を返す
LET rc=0 !完了
FOR i=R TO 1 STEP -1 !N進法R桁の数として
   IF A(i)<N THEN !1～Nで更新する
      LET A(i)=A(i)+1
      FOR j=i+1 TO R !A(i)=A(i+1)= … =A(R)　最初の並び
         LET A(j)=1
      NEXT j
      LET rc=1 !未了
      EXIT SUB
   END IF
NEXT i
END SUB



!最小完全ハッシュ関数 combinadic

!●異なるn個のものから重複を許さずr個を取り出す組合せ

!comb(n,r)通りの組合せパターン ⇔ 0 ～ comb(n,r)-1 の番号

EXTERNAL FUNCTION Comb2Num(A(),N,R) !組合せパターンに番号を付ける　※辞書式順序
LET v=COMB(N,R)-1
FOR j=1 TO R !組合せをビット位置とする
   LET t=N-A(j)
   LET v=v-COMB(t,R-j+1)
NEXT j
LET Comb2Num=v
END FUNCTION

EXTERNAL SUB Num2Comb(h, A(),N,R) !番号から組合せパターンを生成する　※辞書式順序
LET v=COMB(N,R)-h
FOR j=R TO 1 STEP-1 !組合せをビット位置とする
   FOR i=N-1 TO 0 STEP -1 !組合せをビット位置とする
      LET t=COMB(i,j)
      IF v>t THEN EXIT FOR
   NEXT i
   LET v=v-t
   LET A(R-j+1)=N-i !ビット位置(N-i-1)を１とする
NEXT j
END SUB


EXTERNAL FUNCTION Comb2Num2(A(),N,R) !組合せパターンに番号を付ける　※辞書式順序ではない
LET v=0
FOR j=1 TO R !組合せをビット位置とする
   LET t=A(j)-1
   LET v=v+COMB(t,j)
NEXT j
LET Comb2Num2=v
END FUNCTION

EXTERNAL SUB Num2Comb2(h, A(),N,R) !番号から組合せパターンを生成する　※辞書式順序ではない
LET v=h
LET j=R
FOR i=N-1 TO 0 STEP -1 !組合せをビット位置とする
   LET t=COMB(i,j)
   IF t<=v THEN
      LET A(j)=i+1 !ビット位置iを１とする
      LET j=j-1
      LET v=v-t
   END IF
NEXT i
END SUB

EXTERNAL SUB PrevComb(A(),N,R, rc) !辞書式順序で前の組合せを返す
LET rc=0 !完了
FOR i=R TO 2 STEP -1 !並びの差が２以上の位置を探す
   IF A(i)-A(i-1)>1 THEN EXIT FOR
NEXT i
IF i>1 OR (i=1 AND A(1)>1) THEN
   LET A(i)=A(i)-1
   FOR j=i+1 TO R !A(2)=N-R+2< … A(R-1)=N-1<A(R)=N　最後の並び
      LET A(j)=N-R+j
   NEXT j
   LET rc=1 !未了
END IF
END SUB

EXTERNAL SUB NextComb(A(),N,R, rc) !辞書式順序で次の組合せを返す
LET rc=0 !完了
FOR i=R TO 1 STEP -1
   IF A(i)<N-R+i THEN !i～N-R+iで更新する
      LET A(i)=A(i)+1
      FOR j=i+1 TO R !A(i)<A(i+1)< … <A(R)　最初の並び
         LET A(j)=A(j-1)+1
      NEXT j
      LET rc=1 !未了
      EXIT SUB
   END IF
NEXT i
END SUB


!●異なるn個のものから重複を許してr個を取り出す組合せ（重複組合せ）

!comb(n+r-1,r)通りの組合せパターン ⇔ 0 ～ comb(n+r-1,r)-1 の番号

!※Homogeneous Product、nHr=comb(n+r-1,r)=comb(n+r-1,n-1)

EXTERNAL FUNCTION ReptComb2Num(A(),N,R) !組合せパターンに番号を付ける　※辞書式順序
LET v=COMB(N+R-1,R)-1
FOR j=R-1 TO 0 STEP -1
   LET t=N-A(R-j) !0～N-1
   LET v=v-COMB(t+j,j+1) !tHj+1
NEXT j
LET ReptComb2Num=v
!!!別解
!!!FOR j=1 TO R
!!!   LET t=N-A(j) !0～N-1
!!!   LET v=v-COMB(t+R-j,R-j+1) !tHr-j+1
!!!NEXT j
!!!LET ReptComb2Num=v
END FUNCTION

EXTERNAL SUB Num2ReptComb(h, A(),N,R) !番号から組合せパターンを生成する　※辞書式順序
LET v=COMB(N+R-1,R)-h
FOR j=R TO 1 STEP-1 !組合せをビット位置とする
   FOR i=N-1 TO 0 STEP -1 !組合せをビット位置とする
      LET t=COMB(i+j-1,j)
      IF v>t THEN EXIT FOR
   NEXT i
   LET v=v-t
   LET A(R-j+1)=N-i !ビット位置(N-i-1)を１とする
NEXT j
END SUB


EXTERNAL SUB PrevReptComb(A(),N,R, rc) !辞書式順序で前の組合せを返す
LET rc=0 !完了
FOR i=R TO 2 STEP -1 !並びの差が１以上の位置を探す
   IF A(i)-A(i-1)>0 THEN EXIT FOR
NEXT i
IF i>1 OR (i=1 AND A(1)>1) THEN
   LET A(i)=A(i)-1
   FOR j=i+1 TO R !A(i+1)= … A(R-1)=A(R)=N　最後の並び
      LET A(j)=N
   NEXT j
   LET rc=1 !未了
END IF
END SUB

EXTERNAL SUB NextReptComb(A(),N,R, rc) !辞書式順序で次の組合せを返す
LET rc=0 !完了
FOR i=R TO 1 STEP -1
   IF A(i)<N THEN !i～Nで更新する
      LET A(i)=A(i)+1
      FOR j=i+1 TO R !A(i)=A(i+1)= … =A(R)　最初の並び
         LET A(j)=A(j-1)
      NEXT j
      LET rc=1 !未了
      EXIT SUB
   END IF
NEXT i
END SUB

DEF f(x1,x2)=2*x1^2-x1*x2+x2^2-2*x1-3*x2 !非線形関数の最小値
DEF df1(x1,x2)=4*x1-x2-2 !∂f/∂x1
DEF df2(x1,x2)=-x1+2*x2-3 !∂f/∂x2

LET x1=0 !初期値(x1,x2)=(0,0)
LET x2=0

LET cEPS=1E-6 !精度

LET iter=300 !繰り返し回数
FOR i=1 TO iter
   PRINT "x1=";x1, "x2=";x2, "f=";f(x1,x2)

   LET dx1=-df1(x1,x2) !勾配（gradient）
   LET dx2=-df2(x1,x2)
   LET norm2=dx1*dx1+dx2*dx2
   IF norm2<cEPS*cEPS THEN EXIT FOR !勾配が０なら、終了！

   !!!LET alpha=0.1 !ステップ幅（微小な定数による）
   LET alpha=GoldenSelection(x1,x2,dx1,dx2) !直線探索でステップ幅を計算する
   PRINT "α=";alpha !debug

   LET x1=x1+alpha*dx1 !近似解
   LET x2=x2+alpha*dx2
NEXT i

PRINT "繰り返し回数=";i


!１変数関数の最小化
!　目的関数Φk(α)=f(x[k]+α*d[k])を最小にするα∈[0,1]をみつける

FUNCTION GoldenSelection(x1,x2,dx1,dx2) !黄金分割法 GoldenSelection(x,d)
   LET t=(SQR(5)-1)/2 !1/τ、τはτ^2-τ-1=0の解
   LET a=0 !α∈[0,1]
   LET b=1
   LET w=b-a
   LET p=b-t*w !a<p<q< b
   LET q=a+t*w
   LET fxp=f(x1+p*dx1,x2+p*dx2)
   LET fxq=f(x1+q*dx1,x2+q*dx2)
   DO UNTIL w<cEPS !区間が十分小さくなれば終了
      IF fxp>=fxq THEN !区間を[p,b]に縮小する
         LET a=p
         LET w=b-a
         LET p=q
         LET fxp=fxq
         LET q=a+t*w
         LET fxq=f(x1+q*dx1,x2+q*dx2)
      ELSE ![a,q]
         LET b=q
         LET w=b-a
         LET q=p
         LET fxq=fxp
         LET p=b-t*w
         LET fxp=f(x1+p*dx1,x2+p*dx2)
      END IF
   LOOP
   LET GoldenSelection=(a+b)/2 !αの求解
END FUNCTION

END

DEF f(x1,x2)=2*x1^2-x1*x2+x2^2-2*x1-3*x2 !非線形関数の最小値
DEF df1(x1,x2)=4*x1-x2-2 !∂f/∂x1
DEF df2(x1,x2)=-x1+2*x2-3 !∂f/∂x2

LET x1=0 !初期値(x1,x2)=(0,0)
LET x2=0

LET cEPS=1E-6 !精度

LET iter=300 !繰り返し回数
FOR i=1 TO iter
   PRINT "x1=";x1, "x2=";x2, "f=";f(x1,x2)

   LET dx1=-df1(x1,x2) !勾配（gradient）
   LET dx2=-df2(x1,x2)
   LET norm2=dx1*dx1+dx2*dx2
   IF norm2<cEPS*cEPS THEN EXIT FOR !勾配が０なら、終了！

   LET fx=f(x1,x2) !Armijoの基準でステップ幅を計算する
   FOR m=0 TO 10
      LET alpha=1/2^m
      IF f(x1+alpha*dx1,x2+alpha*dx2)<=fx-0.5*alpha*norm2 THEN EXIT FOR
   NEXT m
   PRINT m;alpha !debug

   LET x1=x1+alpha*dx1 !近似解
   LET x2=x2+alpha*dx2
NEXT i

PRINT "繰り返し回数=";i

END

!ジョンソン法 - 順序付け問題

!先に機械M1、次に機械M2を使って生産する製品A,B,C,…の工程順番を決める

LET N=5 !製品の数

DIM T(2,N) !加工時間
!    A,B,C,D,E
DATA 3,5,1,6,7 !M1
DATA 6,2,2,6,5 !M2
MAT READ T

DIM W(2,N)
MAT W=T !copy it

LET cMAX=9999

DIM A(N) !工程の順番
LET Top=1
LET Btm=N
FOR k=1 TO N

   LET s=cMAX !最小のものを探す
   FOR i=1 TO 2
      FOR j=1 TO N
         IF W(i,j)< s THEN
            LET s=W(i,j)
            LET y=i
            LET x=j
         END IF
      NEXT j
   NEXT i
   !!!PRINT y;x;s !debug

   IF y=1 THEN !M1側なら
      LET A(Top)=x !最初へ
      LET Top=Top+1
   ELSE
      LET A(Btm)=x !最後へ
      LET Btm=Btm-1
   END IF

   LET W(1,x)=cMAX !削除して、次へ
   LET W(2,x)=cMAX
   !!!MAT PRINT W; !debug

NEXT k
MAT PRINT A; !debug


DIM M1(N),M2(N) !工程の終了時間
LET s1=0
LET s2=0
FOR i=1 TO N
   LET s1=s1+T(1,A(i)) !つめる
   LET M1(i)=s1

   IF s1>s2 THEN !M1が終了しているかどうかを確認する
      LET s2=s1+T(2,A(i)) !M1の終了を待つ
   ELSE
      LET s2=s2+T(2,A(i)) !つめる
   END IF
   LET M2(i)=s2
NEXT i
MAT PRINT M1;
MAT PRINT M2;

PRINT "最小時間=";MAX(M1(N),M2(N))


!ガントチャートを描く

SET bitmap SIZE 25*20,11*20 !作画領域
SET WINDOW -1,25,-1,11
DRAW grid
SET TEXT JUSTIFY "center","half"
FOR i=1 TO N
   LET u=A(i) !工程

   LET v=T(1,A(i)) !M1側
   LET x2=M1(i)
   LET x1=x2-v
   CALL box(x1,x2,6,8,u,v)

   LET v=T(2,A(i)) !M2側
   LET x2=M2(i)
   LET x1=x2-v
   CALL box(x1,x2,2,4,u,v)
NEXT i
SUB box(x1,x2,y1,y2,u,v) !チャートを描く
   PLOT LINES: x1,y1; x2,y1; x2,y2; x1,y2; x1,y1 !長方形
   PLOT TEXT ,AT (x1+x2)/2,(y1+y2)/2: mid$("ABCDE",u,1)&"("&STR$(v)&")" !工程と時間
END SUB

END

0 回目
 1  2  3  4  5  6  7  8  9  10

 1 回目
 6  5  7  4  8  3  9  2  10  1

 2 回目
 3  8  9  4  2  7  10  5  1  6

 3 回目
 7  2  10  4  5  9  1  8  6  3

 4 回目
 9  5  1  4  8  10  6  2  3  7

 5 回目
 10  8  6  4  2  1  3  5  7  9

 6 回目
 1  2  3  4  5  6  7  8  9  10

!「シャッフル」のシミュレーション

LET N=10 !カードの枚数

DIM shuffle(N) !シャッフルに相当する置換
FOR i=1 TO N
   LET shuffle(i)=INT((N-i*(-1)^i-MOD(N,2))/2)+1
   !LET shuffle(i)=MOD(2*i,N+1) !in shuffle　※Nは２の倍数
   !IF i=N THEN LET shuffle(i)=N ELSE LET shuffle(i)=MOD(2*(i-1),N-1)+1 !out shuffle（1,Nは不変）　※Nは２の倍数
NEXT i
CALL PermPrintOut(shuffle)
CALL PermCyclicPrintOut(shuffle)
!!!STOP


DIM C(N) !カードの束
CALL PermIdentity(C) !カードを整列させる
PRINT "0 回目"
MAT PRINT C;


DIM T(N)
LET cMax=500 !最大回数　※調整が必要
FOR x=1 TO cMax

   CALL PermMultiply(C,shuffle, C) !シャッフル
   PRINT x;"回目"
   MAT PRINT C;

   IF PermIsIdentity(C)<>0 THEN EXIT FOR !元に戻ったら、終了！

NEXT x
IF x>cMax THEN PRINT cMax;"回では元に戻りません。"


END


EXTERNAL SUB PermIdentity(A()) !恒等置換
FOR i=1 TO UBOUND(A)
   LET A(i)=i
NEXT i
END SUB

EXTERNAL FUNCTION PermIsIdentity(A()) !恒等置換か確認する
LET PermIsIdentity=0 !false
FOR i=1 TO UBOUND(A)
   IF A(i)<>i THEN EXIT FUNCTION
NEXT i
LET PermIsIdentity=-1 !true
END FUNCTION

EXTERNAL SUB PermMultiply(A(),B(), AB()) !積AB
DIM T(UBOUND(A))
FOR i=1 TO UBOUND(A)
   LET T(i)=A(B(i)) !※合成写像(AB)(i)=A(B(i))
NEXT i
MAT AB=T
END SUB


!補助ルーチン

EXTERNAL SUB PermPrintOut(A()) !表示する　※標準形（２行ｎ列の行列表記する）
PRINT "┌";
FOR i=1 TO UBOUND(A)
   PRINT USING "###": i;
NEXT i
PRINT " ┐"

PRINT "└";
FOR i=1 TO UBOUND(A)
   PRINT USING "###": A(i);
NEXT i
PRINT " ┘"
END SUB

EXTERNAL SUB PermCyclicPrintOut(A()) !巡回置換の積で表示する　※(,,…,)(,,…,)…(,,…,)
DIM T(UBOUND(A)) !作業用
CALL PermIdentity(T)

FOR i=1 TO UBOUND(A)
   IF T(i)>0 THEN !含まれないなら
      LET L=0 !長さ
      PRINT "(";

      LET k=i !最初の値
      DO !リストを作成する　※(4 4)などの長さ１の置換も含む
         PRINT k;
         LET L=L+1

         LET T(k)=-1 !巡回に含まれる

         LET k=A(k) !次をトレースする
      LOOP UNTIL k=i !一巡するまで

      PRINT ")"; !次の組へ
      PRINT " 長さ=";L
   END IF
NEXT i
PRINT
END SUB