% ----------------------------------------------------------- % % Basics of Matrix and Linear Algebra by Prolog % ----------------------------------------------------------- % % file: matrix1b.pl % imported: (math1.pl version: 21 Feb 2003) % edited: 25-26 Feb 2003.(matrix1.pl) % modified: 18 Apr 2003. case_of_diagonal(4,Pxx,_,_,_,_). % modified: 8-9 Jun 2003. symbolic equations. % modified: 20 Aug 2003. matrix_is with modulo. % modified: 20 Aug 2003. the program has rewrited with % some new predicates len/2, kth/3, and rev/2 each of which % is the alternative of SWI system predicates length/2, % nth1/3, and reverse/2 respectively. :- dynamic matrix /2. :- dynamic cell /3. % ----------------------------------------------------------- % % definitions of matrix representation % ----------------------------------------------------------- % % cell(1,1,1). cell(1,2,2). cell(2,1,3). cell(2,2,4). matrix(i(2),[[1,0],[0,1]]). matrix(i(3),[[1,0,0],[0,1,0],[0,0,1]]). matrix(ma,[[0,-1],[1,0]]). matrix(mb,[[1,-1],[2,0]]). matrix(mc,[[3,-6,-2],[1,-2,-1],[-2,6,3]]). matrix(p,[[1,3,2],[0,1,1],[1,0,-2]]). matrix(q,[[2,-6,-1],[-1,4,1],[1,-3,-1]]). matrix(r,[[1,-1,-1],[-1,1,1],[1,-1,-1]]). matrix(s,[[1,1,1],[2,1,-1],[3,-1,1]]). matrix(t,[[1,1,1],[1,1,1],[3,-1,1]]). matrix(n2sq,[[a,b],[c,d]]). matrix(x,[[x,z],[y,w]]). % identity matrix matrix(i(N),I):- integer(N), N > 3, matrix_of_size(I,N,N), findall(R, ( kth(K,I,R), findall(C, ( kth(L,R,C), (L = K -> C = 1; C = 0) ), R) ), I),!. % scalar multiplication % modified: 8 June 2003. symbolic case. matrix(C * X,B):- %eval_number(C,_), matrix(X,A), matrix_of_size(A,N,M), matrix_of_size(B,N,M), findall(Br, ( kth(K,A,Ar), kth(K,B,Br), findall(Brc, ( kth(L,Ar,Arc), kth(L,Br,Brc), evaluate_if_numerical(Brc is C * Arc) ), Br) ), B), !. % zero matrix matrix(o(N,M),O):- matrix(z(N,M,0),O). % uniform bias matrix(z(N,M,C),O):- \+ C == cell, matrix_of_size(O,N,M), findall(Ok, ( kth(_K,O,Ok), findall(C, ( kth(_L,Ok,C) ), Ok) ), O),!. % indirect definition via cells. matrix(z(N,M,cell),O):- X =cell(_,_,_), clause(X,_), matrix_of_size(O,N,M), findall(Ok, ( kth(K,O,Ok), findall(C, ( kth(L,Ok,C), cell(K,L,C) ), Ok) ), O),!. % numerical evaluation : a matrix version of `is /2' % modified: 8 June 2003. for cases of symbolic matrix_is(C, A):- matrix(C is A). matrix(C is X):- \+ var(X), (matrix(X,A);X=A), matrix_of_size(A,N,M), matrix_of_size(C,N,M), findall(Ck, ( kth(K,A,Ak), kth(K,C,Ck), findall(Ckl, ( kth(L,Ak,Akl), evaluate_if_numerical(Ckl is Akl), kth(L,Ck,Ckl) ), Ck) ), C),!. % mod % added: 20 Aug 2003. in order for design0.pl matrix_mod(C is (A mod 2^2)):- !, C = A. matrix_mod(C is (A mod P)):- matrix_mod(C, A, mod(P)). matrix_mod(C, A, mod(P)):- \+ var(A), integer(P), (matrix(X,A);X=A), matrix_of_size(A,N,M), matrix_of_size(C,N,M), findall(Ck, ( kth(K,A,Ak), kth(K,C,Ck), findall(Ckl, ( kth(L,Ak,Akl), Ckl is Akl mod P, kth(L,Ck,Ckl) ), Ck) ), C),!. % ----------------------------------------------------------- % % basic calculations of matrix and linear algebra % ----------------------------------------------------------- % % modified: 8 June 2003. % matrix_add(X + Y = Z):- (matrix(X,A);X=A), (matrix(Y,B);Y=B), matrix_of_size(A,N,M), matrix_of_size(B,N,M), findall(ZK, ( kth(K,A,AK), kth(K,B,BK), findall(ZKL, ( kth(L,AK,AKL), kth(L,BK,BKL), evaluate_if_numerical(AKL1 is AKL), evaluate_if_numerical(BKL1 is BKL), evaluate_if_numerical(ZKL is AKL1 + BKL1) ), ZK) ), Z). matrix_negative(- Y = Z):- matrix(Y,B), matrix_of_size(B,_N,_M), findall(ZK, ( kth(_K,B,BK), findall(ZKL, ( kth(_L,BK,BKL), evaluate_if_numerical(ZKL is (-1)*BKL) ), ZK) ), Z). matrix_subtract(X - Y = Z):- matrix_negative(- Y = Z1), retractall(matrix(temp,_)), assert(matrix(temp,Z1)), matrix_add(X + temp = Z). % matrix generator % ----------------------------------------------------------- % matrix_of_size(A,N,M):- len(A,N), kth(1,A,A1), len(A1,M), findall(B, ( member(B,A), len(B,M) ), A). % restriction by rows and columns. % ----------------------------------------------------------- % restricted_matrix(X,[B],[C],[[D]]):- (matrix(X,A);X=A), matrix_of_size(A,_N0,_M0), kth(B,A,Ab), kth(C,Ab,D). restricted_matrix(X,B,C,D):- % X: original matrix % B: restriction for rows % C: restriction for columns % D: restricted matrix B \= [_], C \= [_], (matrix(X,A);X=A), matrix_of_size(A,N0,M0), len(B,N), len(C,M), N =< N0, M =< M0, matrix_of_size(D,N,M), findall(DK, ( member(K,B), kth(K,A,AK), findall(AKL, ( member(L,C), kth(L,AK,AKL) ), DK) ), D). % vertical / horizontal projections % ----------------------------------------------------------- % matrix_rows(A,[R],[B]):- kth(R,A,B). matrix_rows(A1,R,B):- \+ R = [_], (matrix(A1,A);A=A1), matrix_of_size(A,_N,M), anum_seq1(Q,M), restricted_matrix(A,R,Q,B). matrix_cols(A1,C,B):- (matrix(A1,A);A=A1), matrix_of_size(A,N,_M), anum_seq1(Q,N), restricted_matrix(A,Q,C,B). % index for matrix matrix_index(X,A,B,[C,D]):- (matrix(X,A);X=A), matrix_of_size(A,_N0,_M0), restricted_matrix(A,[C],[D],[[B]]). % cf., collect the L(i)-th element from each row % acording to a choice list L. matrix_collect(X,A,B,L):- matrix(X,A), index_of_tuple(A,B,L). % transposition % ----------------------------------------------------------- % matrix_transpose(Y -> Z):- (matrix(Y,A);Y=A), matrix_of_size(A,N,M), matrix_of_size(Z,M,N), findall(ZK, ( kth(K,Z,ZK), findall(ZKL, ( kth(L,ZK,ZKL), matrix_index(Y,A,ZKL,[L,K]) ), ZK) ), Z). % multiplication % ----------------------------------------------------------- % % modified: 8 June 2003. matrix_multiply(X * Y = Z):- \+ number(X), (matrix(X,A);X=A), (matrix(Y,B);Y=B), matrix_of_size(A,N,L0), matrix_of_size(B,L0,M), matrix_of_size(Z,N,M), matrix_transpose(Y -> C),!, findall(ZK, ( kth(K,Z,ZK), findall(ZKL, ( kth(L,ZK,ZKL), matrix_rows(A,[K],[AK]), matrix_rows(C,[L],[CL]), eq_sum_product(AK,CL,_,ZKL,_) ), ZK) ), Z). matrix_multiply(X * Y = Z):- eval_number(X,_), (matrix(Y,B);Y=B), matrix_of_size(B,_N,M), matrix(X * i(M),Ix), matrix_multiply(Ix * Y = Z1), matrix(Z is Z1). matrix_multiply_series([X],X,1):- matrix(X,A);matrix(A,X). matrix_multiply_series([X|Y],Z,K):- matrix_multiply_series(Y,Z1,K1), K is K1 + 1, matrix_multiply(X * Z1 = Z). matrix_powered(X ^ Y = Z):- integer(Y), matrix(X,_A), bag0(B,[X],Y), matrix_multiply_series(B,Z,Y). % expanded matrix % ----------------------------------------------------------- % matrix_expand(X | Y = C):- (matrix(X,A);X=A), (matrix(Y,B);Y=B), matrix_of_size(A,N,M), matrix_of_size(B,N,M1), M2 is M + M1, matrix_of_size(C,N,M2), findall(CK, ( kth(K,A,AK), kth(K,B,BK), append(AK,BK,CK) ), C), !. % modified: 8 Jun 2003. % ----------------------------------------------------------- % row_vector_multiply(C*B=D):- len(C,M), len(B,M), eq_sum_product(B,C,D,_,_). row_vector_multiply(C0*B=D):- evaluate_if_numerical(C is C0), len(B,M), bag0(CL,[C],M), eq_sum_product(B,CL,D,_,_). % % the inverse of regular matrix. % ----------------------------------------------------------- % matrix_inverse(P ^ -1 = Q):- (matrix(P,B);P=B), matrix_of_size(B,N,M), len(Xs,N), eliminate_3(P,Xs,_,V), %matrix_of_size(Q,N,M), %matrix_expand(i(N)|Q = V). % good sense. but... anum_seq1(Ys1,N), matrix(z(1,M,M),O), matrix_add([Ys1]+O=R), matrix([Ys] is R), matrix_cols(V,Ys,Q), !. /* % a sample execution. ?- matrix_inverse(p^(-1)=Q),matrix_multiply(p*Q=E),matrix(F is E). Q = [[2, -6, -1], [-1, 4, 1], [1, -3, -1]] E = [[2*1+3* -1+1*2, 2* -3+3*4+1* -6, 2* -1+3*1+1* -1], [1*1+1* -1+0*2, 1* -3+1*4+0* -6, 1* -1+1*1+0* -1], [-2*1+0* -1+1*2, -2* -3+0*4+1* -6, -2* -1+0*1+1* -1]] F = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] */ % ----------------------------------------------------------- % % Gauss elimination % ----------------------------------------------------------- % % % row operations % ----------------------------------------------------------- % % modified: 8 Jun 2003 row_operation(P,0,[_,_], _, P). row_operation(P,C0,[K,L], A + C * B, FF):- % +P: a matrix % K,A: target row (K-th row) % L,B: pivot row (L-th row) % +C: constant % E: outcome evaluate_if_numerical(C is C0), \+ C = 0, (matrix(P,Q);Q=P), matrix_of_size(Q,N,M), matrix_of_size(E,N,M), kth(K,Q,A), kth(L,Q,B), %matrix_rows(Q,[K],[A]), %matrix_rows(Q,[L],[B]), %matrix(A1 is A), %matrix(B1 is B), row_vector_multiply(C*B=CB), matrix_add([A]+[CB]=XXX), matrix([ACB] is XXX), replace(K/N,Q,ACB,E), nl,write(row_operation(row(K)=row(K)+const(C)*row(L))), matrix(FF is E), nl,write((P)), nl,write((->)), nl,write((FF)). row_operation_series(1,P,[[C,K,L]],[A + C * B],[E]):- row_operation(P,C,[K,L], A + C * B, E). row_operation_series(K,P,[[C,K,L]|X], [A + C * B|Y], [E|[F|G]]):- row_operation_series(K1,P,X,Y,[F|G]), K is K1 + 1, row_operation(F,C,[K,L], A + C * B, E). expand_to_eliminate(A,B):- (matrix(A,P);P=A), matrix_of_size(P,N,N), matrix(i(N),I), matrix_expand(P|I=B). % % to make the diagonal element c(X,X)=1. % ----------------------------------------------------------- % % level 2 diagonalize_2(P,Q,[],[],[]):- (matrix(P,Q);P=Q), matrix_of_size(Q,_N,_M). %expand_to_eliminate(P,B). diagonalize_2(P,Q,[E|F],[G|H],[X|Y]):- diagonalize_2(P,E,F,H,Y), len(E,N), (len(Y,N)->!,fail;true), matrix_index(E,_,Pxx,[X,X]), \+ member(X,Y), diagonalize_1((E,X,Pxx),G,Q,_). % base level diagonalize_1((B,X,Pxx),(L,C,A),P1,E):- % X,Pxx: target of diagnolization matrix_index(B,B,Pxx,[X,X]), case_of_diagonal(Q,Pxx,X,B,C,L), row_operation(B,C,[X,L], A, P1), %~~~~~~~~~~~~~~~~~~~~~~~! matrix(E is P1), nl,write(case_dialog_1((Q,X,X))). % modified: 8 Jun 2003. case_of_diagonal(0,Pxx,X,_,C,X):- ( (Pxx=..[_|Sxx],Sxx\=[],member(X,Sxx),\+ number(X)); (Pxx=..[Sxx],\+ number(Sxx)) ), C1 = 1/Pxx-1, evaluate_if_numerical(C is C1), !. case_of_diagonal(1,Pxx,X,_,0,X):- Vxx is Pxx, Vxx is 1, !. case_of_diagonal(2,Pxx,X,_,-2,X):- Vxx is Pxx, Vxx is -1, !. case_of_diagonal(3,Pxx,X,_,C,X):- Vxx is Pxx, (Vxx > 0; Vxx < 0), C is (1-Pxx)/Pxx, !. case_of_diagonal(4,Pxx,X,B,1/Plx,L):- Vxx is Pxx, Vxx is 0, matrix_index(B,B,Plx,[L,X]), L > X, % added: 18 Apr 2003. Vlx is Plx, (Vlx > 0; Vlx < 0), !. % % elimination of non diagonal elements. % ----------------------------------------------------------- % eliminate_3(P,[X],[Q],V):- expand_to_eliminate(P,PE), len(PE,N), matrix_index(PE,PE,_Pxx,[X,X]), len(Xs,N), eliminate_2(PE,PE,X,Xs,_,Q), matrix(V is Q). eliminate_3(P,[X|Y],[Q|R],V):- eliminate_3(P,Y,[Q1|R],_V1), len(Q1,N), (len(Y,N)->!,fail;true), matrix_index(Q1,Q1,_Pxx,[X,X]), \+ member(X,Y), len(Xs,N), eliminate_2(Q1,Q1,X,Xs,_,Q), matrix(V is Q). % % eliminate non-zero in the X-th column. % eliminate_2(P,B,X,[K|G],[Q|H],V):- eliminate_2(P,B,X,G,[Q1|H],_V1), len(Q1,N), (len(G,N)->!,fail;true), matrix_index(Q1,Q1,Pkx,[K,X]), \+ member(K,G), case_of_elimination(_EC,Pkx,C), row_operation(Q1,C,[K,X], _A, Q), matrix(V is Q), nl,write(elim2(case_of_elim(_EC),[row:K,col:X])). % % start with diagonals of 1 (if numerical). % modified: 8 Jun 2003. eliminate_2(P,B,X,[X],[Q],V):- eliminate_1(P,B,Q,V), matrix_index(Q,Q,Pxx,[X,X]), case_of_diagonal(Case,Pxx,X,_,C0,X), member(Case,[0,1]), nl,write(elim2(case_of_diag(1),[row:X,col:X],C0)). eliminate_1(P,B,Q,V):- (matrix(P,B);P=B), matrix_of_size(B,N,_M), anum_seq1(Y,N), diagonalize_2(P,Q,_,_,Y), matrix(V is Q). % added: 8 Jun 2003. case_of_elimination(0,Pkx,Qkx):- ( (Pkx=..[_|Skx],Skx\=[],member(X,Skx),\+ number(X)); (Pkx=..[Skx],\+ number(Skx)) ), evaluate_if_numerical(Qkx is (-Pkx)), !. case_of_elimination(1,Pkx,-1):- Vkx is Pkx, Vkx is 1, !. case_of_elimination(2,Pkx,1):- Vkx is Pkx, Vkx is -1, !. case_of_elimination(3,Pkx,Qkx):- Vkx is Pkx, (Vkx > 0; Vkx < 0), Qkx is - Vkx, !. case_of_elimination(4,Pkx,0):- Vkx is Pkx, Vkx is 0, !. % ----------------------------------------------------------- % % Arithmetic and so on including probabilistic operators % ----------------------------------------------------------- % % edited: 26 Feb 2003. (imported from: math1.pl) % % evaluate the numeric value 'if' it is numeric equation % ----------------------------------------------------------- % % added: 8-9 June 2003 evaluate_if_numerical(X1 is X):- functor(X,Functor,Arity), ( (member(Functor,[+,-,*,/]),Arity=2); (member(Functor,[-]),Arity=1) ), X=..[Functor|Args], ( forall( member(Y,Args), number(Y)) -> X1 is X ; ( eval_equation(step1,X,XX0), eval_equation(step1b,XX0,XX1), eval_equation(step2,XX1,XX2), eval_equation(step3,XX2,X1) ) ), !. evaluate_if_numerical(X is X):- functor(X,Functor,Arity), \+ member(Functor,[+,-,*,/]), Arity=0, !. % % evaluation of a nummerical value % ----------------------------------------------------------- % %eval_equation(step1,X,X1) eval_equation(step1,_/0, infinite). eval_equation(step1,0/_, 0). eval_equation(step1,A/A, 1). eval_equation(step1,(A/B)*(B/A), 1). eval_equation(step1,A/(C/A), 1/C). %------- eval_equation(step1,0*_Arg2, 0). eval_equation(step1,_Arg1*0, 0). eval_equation(step1,0+Arg2, Arg2). eval_equation(step1,0-Arg2, Arg2). eval_equation(step1,Arg1+0, Arg1). eval_equation(step1,Arg1-0, Arg1). eval_equation(step1,A+(C-A), C). eval_equation(step1,A-(C+A), C). eval_equation(step1,(-A)+(A+C), C). eval_equation(step1,(-A)+(C+A), C). eval_equation(step1,(B+C)-B, C). eval_equation(step1,(C+B)-B, C). eval_equation(step1,(C-B)+B, C). eval_equation(step1,((-B)+C)+B, C). %------- eval_equation(step1,1*B, B). eval_equation(step1,B*1, B). eval_equation(step1,(-1)*B, -B). eval_equation(step1,B*(-1), -B). eval_equation(step1,(-A)+A, 0). eval_equation(step1,A-A, 0). eval_equation(step1,(-A)/A, -1). eval_equation(step1,A/(-A), -1). %------- eval_equation(step1,(-1/A)*A, -1). eval_equation(step1,A*(-1/A), -1). eval_equation(step1,-(1/A)*A, -1). eval_equation(step1,A*(-(1/A)), -1). eval_equation(step1,((-1)/A)*A, -1). eval_equation(step1,A*((-1)/A), -1). eval_equation(step1,(1/A)*(-A), -1). eval_equation(step1,(-A)*(1/A), -1). eval_equation(step1,(1/A)*A, 1). eval_equation(step1,A*(1/A), 1). eval_equation(step1,(-C/A)*A, -C). eval_equation(step1,A*(-C/A), -C). eval_equation(step1,(-(C/A))*A, -C). eval_equation(step1,A*(-(C/A)), -C). eval_equation(step1,((-C)/A)*A, -C). eval_equation(step1,A*((-C)/A), -C). eval_equation(step1,(C/A)*(-A), -C). eval_equation(step1,(-A)*(C/A), -C). eval_equation(step1,(C/A)*A, C). eval_equation(step1,A*(C/A), C). %------- eval_equation(step1,(-C/A)*(B/D), -(B*C)/(A*D)). eval_equation(step1,(B/D)*(-C/A), -(B*C)/(A*D)). eval_equation(step1,(-(C/A)*(B/D)), -(B*C)/(A*D)). eval_equation(step1,(-(B/D)*(C/A)), -(B*C)/(A*D)). eval_equation(step1,(-(C/A))*(B/D), -(B*C)/(A*D)). eval_equation(step1,(B/D)*(-(C/A)), -(B*C)/(A*D)). eval_equation(step1,((-C)/A)*(B/D), -(B*C)/(A*D)). eval_equation(step1,(B/D)*((-C)/A), -(B*C)/(A*D)). eval_equation(step1,(C/A)*(-B/D), -(C*B)/(A*D)). eval_equation(step1,(-B/D)*(C/A), -(C*B)/(A*D)). eval_equation(step1,(C/A)*((-B)/D), -(C*B)/(A*D)). eval_equation(step1,((-B)/D)*(C/A), -(C*B)/(A*D)). eval_equation(step1,(C/A)*(A/D), C/D). eval_equation(step1,(A/D)*(C/A), C/D). eval_equation(step1,(C/A)*(B/D), (C*B)/(A*D)). %------- eval_equation(step1,(-C/A)*(A/D), -C/D). eval_equation(step1,(A/D)*(-C/A), -C/D). eval_equation(step1,(-(C/A)*(A/D)), -C/D). eval_equation(step1,(-(A/D)*(C/A)), -C/D). eval_equation(step1,(-(C/A))*(A/D), -C/D). eval_equation(step1,(A/D)*(-(C/A)), -C/D). eval_equation(step1,((-C)/A)*(A/D), -C/D). eval_equation(step1,(A/D)*((-C)/A), -C/D). eval_equation(step1,(C/A)*(-A/D), -C/D). eval_equation(step1,(-A/D)*(C/A), -C/D). eval_equation(step1,(C/A)*((-A)/D), -C/D). eval_equation(step1,((-A)/D)*(C/A), -C/D). %------- eval_equation(step1,(-1/A)*B, -B/A). eval_equation(step1,B*(-1/A), -B/A). eval_equation(step1,(-(1/A))*B, -B/A). eval_equation(step1,B*(-(1/A)), -B/A). eval_equation(step1,((-1)/A)*B, -B/A). eval_equation(step1,B*((-1)/A), -B/A). eval_equation(step1,(1/A)*(-B), -B/A). eval_equation(step1,(-B)*(1/A), -B/A). eval_equation(step1,(1/A)*B, B/A). eval_equation(step1,B*(1/A), B/A). eval_equation(step1,(-C/A)*B, -C*B/A). eval_equation(step1,B*(-C/A), -C*B/A). eval_equation(step1,(-(C/A))*B, -C*B/A). eval_equation(step1,B*(-(C/A)), -C*B/A). eval_equation(step1,((-C)/A)*B, -C*B/A). eval_equation(step1,B*((-C)/A), -C*B/A). eval_equation(step1,(C/A)*(-B), -C*B/A). eval_equation(step1,(-B)*(C/A), -C*B/A). eval_equation(step1,(C/A)*B, C*B/A). eval_equation(step1,B*(C/A), C*B/A). %------- eval_equation(step1,Arg1+(-B), Arg1-B). eval_equation(step1,Arg1-(-B), Arg1+B). eval_equation(step1,Arg1+B, Arg1-C):- number(B),C is abs(B),C>B. eval_equation(step1,Arg1-B, Arg1+C):- number(B),C is abs(B),C>B. %------- eval_equation(step1,A/(C*A), 1/C). eval_equation(step1,A/(A*C), 1/C). eval_equation(step1,(A*C)/A, C). eval_equation(step1,(A*C)/A, C). %------- eval_equation(step1,(-A)/(C*A), -1/C). eval_equation(step1,(-A)/(A*C), -1/C). eval_equation(step1,A/(C*(-A)), -1/C). eval_equation(step1,A/((-A)*C), -1/C). eval_equation(step1,(A*C)/(-A), -C). eval_equation(step1,(A*C)/(-A), -C). eval_equation(step1,((-A)*C)/A, -C). eval_equation(step1,((-A)*C)/A, -C). %------- eval_equation(step1,(A*B)/(C/A), B/C). eval_equation(step1,(B*A)/(C/A), B/C). eval_equation(step1,A/(C/(A*B)), B/C). eval_equation(step1,A/(C/(B*A)), B/C). eval_equation(step1,A/(C/B), A*B/C). eval_equation(step1,(A/B)/C, A/(B*C)). eval_equation(step1,A/B/C, A/(B*C)). eval_equation(step1,(-A/B/C), (-A/(B*C))). %------- eval_equation(step1,((-A)*B)/(C/A), -B/C). eval_equation(step1,(B*A)/(C/(-A)), -B/C). eval_equation(step1,(-A)/(C/(A*B)), -B/C). eval_equation(step1,(-A)/(C/(B*A)), -B/C). eval_equation(step1,(A*B)/(C/(-A)), -B/C). eval_equation(step1,(B*(-A))/(C/A), -B/C). eval_equation(step1,A/(C/((-A)*B)), -B/C). eval_equation(step1,A/(C/(B*(-A))), -B/C). %------- eval_equation(step1,A+(C*A), (C+1)*A). eval_equation(step1,A+(A*C), (C+1)*A). eval_equation(step1,(A*C)+A, (C+1)*A). eval_equation(step1,(A*C)+A, (C+1)*A). eval_equation(step1,A-(C*A), -(C-1)*A). eval_equation(step1,A-(A*C), -(C-1)*A). eval_equation(step1,(C*a)-A, (C-1)*A). eval_equation(step1,(A*C)-A, (C-1)*A). %------- eval_equation(step1,(-A)+(C*A), (C-1)*A). eval_equation(step1,(-A)+(A*C), (C-1)*A). eval_equation(step1,(A*C)+(-A), (C-1)*A). eval_equation(step1,(A*C)+(-A), (C-1)*A). eval_equation(step1,(-A)-(C*A), -(C+1)*A). eval_equation(step1,(-A)-(A*C), -(C+1)*A). eval_equation(step1,(A*C)-(-A), (C+1)*A). eval_equation(step1,(A*C)-(-A), (C+1)*A). %------- eval_equation(step1,(-A)*B, -(A*B)). eval_equation(step1,A*(-B), -(A*B)). eval_equation(step1,(-A)/B, -A/B). eval_equation(step1,A/(-B), -A/B). eval_equation(step1,(-A)/(-B), A/B). eval_equation(step1,(-(-A/B)), A/B). %------- eval_equation(step1,X,X). eval_equation(step1b,A*B, C):- eval_equation(step1,A, AX), eval_equation(step1,B, BX), eval_equation(step1,AX*BX, C). eval_equation(step1b,A-B, C):- eval_equation(step1,A, AX), eval_equation(step1,B, BX), eval_equation(step1,AX-BX, C). eval_equation(step1b,A/B, C):- eval_equation(step1,A, AX), eval_equation(step1,B, BX), eval_equation(step1,AX/BX, C). eval_equation(step1b,A+B, C):- eval_equation(step1,A, AX), eval_equation(step1,B, BX), eval_equation(step1,AX+BX, C). eval_equation(step1b,X,X). eval_equation(step2,X,X1):- ( member(X, [ (C*A)+(D*A), (C*A)+(A*D), (A*C)+(D*A), (A*C)+(A*D) ] ) -> (eval_equation(step1,C+D,Eq1),X1 = Eq1*A) ); %------- ( member(X, [ (C*(-A))+(D*A), (C*(-A))+(A*D), ((-A)*C)+(D*A), ((-A)*C)+(A*D) ] ) -> (eval_equation(step1,C-D,Eq1),X1 = Eq1*A) ); %------- ( member(X, [ (C*A)+(D*(-A)), (C*A)+((-A)*D), (A*C)+(D*(-A)), (A*C)+((-A)*D) ] ) -> (eval_equation(step1,C-D,Eq1),X1 = Eq1*A) ); %------- ( member(X, [ (C*A)-(D*A), (C*A)-(A*D), (A*C)-(D*A), (A*C)-(A*D) ] ) -> (eval_equation(step1,C-D,Eq1),X1 = Eq1*A) ); %------- ( member(X, [ (C*(-A))-(D*A), (C*(-A))-(A*D), ((-A)*C)-(D*A), ((-A)*C)-(A*D) ] ) -> (eval_equation(step1,(-(C+D)),Eq1),X1 = Eq1*A) ); %------- ( member(X, [ (C*A)-(D*(-A)), (C*A)-((-A)*D), (A*C)-(D*(-A)), (A*C)-((-A)*D) ] ) -> (eval_equation(step1,(-(C+D)),Eq1),X1 = Eq1*A) ); %------- ( member(X, [ (C*A)/(D*A), (C*A)/(A*D), (A*C)/(D*A), (A*C)/(A*D) ] ) -> (eval_equation(step1,(C/D),Eq1),X1 = Eq1*A) ); %------- ( member(X, [ (C*(-A))/(D*A), (C*(-A))/(A*D), ((-A)*C)/(D*A), ((-A)*C)/(A*D) ] ) -> (eval_equation(step1,(-(C/D)),Eq1),X1 = Eq1*A) ); %------- ( member(X, [ (C*A)/(D*(-A)), (C*A)/((-A)*D), (A*C)/(D*(-A)), (A*C)/((-A)*D) ] ) -> (eval_equation(step1,(-(C/D)),Eq1),X1 = Eq1*A) ); %------- X1=X. eval_equation(step4,Num/Den,X1):- eval_equation(step1,Num,NumEV), eval_equation(step1,Den,DenEV), eval_equation(step1,NumEV/DenEV,X1). eval_equation(step5,(Num1+Num2)/Den,X1):- eval_equation(step1,Num1,NumEV1), eval_equation(step1,Num2,NumEV2), eval_equation(step1,NumEV1+NumEV2,NumEV), eval_equation(step1,Den,DenEV), eval_equation(step1,NumEV/DenEV,X1). eval_equation(step3,X,X1):- /* (X=C+(D/A) -> (eval_equation(step1,(C*A),E),eval_equation(step4,[(E+D)/A],X1))); (X=C+(-D/A) -> (eval_equation(step1,(C*A),E),eval_equation(step4,[(E-D)/A],X1))); (X=C+(-(D/A)) -> (eval_equation(step1,(C*A),E),eval_equation(step4,[(E-D)/A],X1))); (X=C+((-D)/A) -> (eval_equation(step1,(C*A),E),eval_equation(step4,[(E-D)/A],X1))); (X=C+(D/(-A)) -> (eval_equation(step1,(C*A),E),eval_equation(step4,[(E-D)/A],X1))); */ %------- (X=(C/A)+(D/A) -> eval_equation(step4,[(C+D)/A],X1)); (X=(-C/A)+(-D/A) -> eval_equation(step4,[(C+D)/A],X1)); (X=(C/A)+(-D/A) -> eval_equation(step4,[(C-D)/A],X1)); (X=(C/A)+(-(D/A)) -> eval_equation(step4,[(C-D)/A],X1)); (X=(C/A)+((-D)/A) -> eval_equation(step4,[(C-D)/A],X1)); (X=(C/A)+(D/(-A)) -> eval_equation(step4,[(C-D)/A],X1)); (X=(-C/A)+(D/A) -> eval_equation(step4,[(D-C)/A],X1)); (X=(-(C/A))+(D/A) -> eval_equation(step4,[(D-C)/A],X1)); (X=(-(C)/A)+(D/A) -> eval_equation(step4,[(D-C)/A],X1)); (X=(C/(-A))+(D/A) -> eval_equation(step4,[(D-C)/A],X1)); %------- (X=(C/A)-(D/A) -> eval_equation(step4,[(C-D)/A],X1)); (X=(-C/A)-(-D/A) -> eval_equation(step4,[(D-C)/A],X1)); (X=(-C/A)-(D/A) -> eval_equation(step4,[(-C-D)/A],X1)); (X=(C/A)-(-D/A) -> eval_equation(step4,[(C-D)/A],X1)); (X=(-(C/A))-(D/A) -> eval_equation(step4,[(-C-D)/A],X1)); (X=((-C)/A)-(D/A) -> eval_equation(step4,[(-C-D)/A],X1)); (X=(C/(-A))-(D/A) -> eval_equation(step4,[(-C-D)/A],X1)); (X=(C/A)-(-(D/A)) -> eval_equation(step4,[(C-D)/A],X1)); (X=(C/A)-((-D)/A) -> eval_equation(step4,[(C-D)/A],X1)); (X=(C/A)-(D/(-A)) -> eval_equation(step4,[(C-D)/A],X1)); %------- (X=(C/A)+(D/B) -> eval_equation(step5,[(C*D+D*A)/AB],X1)); (X=(-C/A)+(D/B) -> eval_equation(step5,[((-C*B)+D*A)/AB],X1)); (X=(-(C/A))+(D/B) -> eval_equation(step5,[((-C*B)+D*A)/AB],X1)); (X=((-C)/A)+(D/B) -> eval_equation(step5,[((-C*B)+D*A)/AB],X1)); (X=(C/(-A))+(D/B) -> eval_equation(step5,[((-C*B)+D*A)/AB],X1)); (X=(C/A)+(-D/B) -> eval_equation(step5,[(C*B-D*A)/AB],X1)); (X=(C/A)+(-(D/B)) -> eval_equation(step5,[(C*B-D*A)/AB],X1)); (X=(C/A)+((-D)/B) -> eval_equation(step5,[(C*B-D*A)/AB],X1)); (X=(C/A)+(D/(-B)) -> eval_equation(step5,[(C*B-D*A)/AB],X1)); %------- (X=(C/A)-(D/B) -> eval_equation(step5,[(C*D-D*A)/AB],X1)); (X=(-C/A)-(D/B) -> eval_equation(step5,[((-C*B)-D*A)/AB],X1)); (X=(-(C/A))-(D/B) -> eval_equation(step5,[((-C*B)-D*A)/AB],X1)); (X=((-C)/A)-(D/B) -> eval_equation(step5,[((-C*B)-D*A)/AB],X1)); (X=(C/(-A))-(D/B) -> eval_equation(step5,[((-C*B)-D*A)/AB],X1)); (X=(C/A)-(-D/B) -> eval_equation(step5,[(C*B+D*A)/AB],X1)); (X=(C/A)-(-(D/B)) -> eval_equation(step5,[(C*B+D*A)/AB],X1)); (X=(C/A)-((-D)/B) -> eval_equation(step5,[(C*B+D*A)/AB],X1)); (X=(C/A)-(D/(-B)) -> eval_equation(step5,[(C*B+D*A)/AB],X1)); %------- X1=X. % % evaluation of a nummerical value % ----------------------------------------------------------- % eval_number(X,X):- number(X),!. % % max,min % ----------------------------------------------------------- % max_of(X,[X]). max_of(Z,[X|Y]):- max_of(Z1,Y), (X > Z1 -> Z=X; Z=Z1). min_of(X,[X]). min_of(Z,[X|Y]):- min_of(Z1,Y), (X < Z1 -> Z=X; Z=Z1). % count frequency of occurence of the specified value of variable, M. % ----------------------------------------------------------- % % note: Both of M and L have to be specified. counter(N,M,L):- len(L,_), findall(M,member(M,L),Mx), len(Mx,N). % sum % ----------------------------------------------------------- % sum([],0). sum([X|Members],Sum):- sum(Members,Sum1), %number(X), Sum is Sum1 + X. % % product % ----------------------------------------------------------- % product([],1). product([X|Members],Z):- product(Members,Z1), %number(X), Z is Z1 * X. % % weighted sum % ----------------------------------------------------------- % product_sum([],[],[],0). product_sum([P|Q],[A|B],[E|F],V):- len(Q,N), len(B,N), product_sum(Q,B,F,V1), E is P * A, V is V1 + E. % % product sum value with equational % ----------------------------------------------------------- % product_sum_eq([],[],[],0,0). product_sum_eq([P|Q],[A|B],[E|F],V,Vq):- len(Q,N), len(B,N), product_sum_eq(Q,B,F,V1,Vq1), Eq = (P) * A, E is Eq, (Vq1=0 -> Vq = Eq; Vq = Vq1 + Eq), V is V1 + E. % % symbolic representation of eq_sum_product (product_sum). % added: 8-9 June 2003 % ----------------------------------------------------------- % eq_sum_product([],[],[],0,0). eq_sum_product([P|Q],[A|B],[E|F],V,Vq):- len(Q,N), len(B,N), eq_sum_product(Q,B,F,V1,Vq1), evaluate_if_numerical(EP is P), evaluate_if_numerical(EA is A), evaluate_if_numerical(E is EP * EA), (Vq1=0 -> Vq = EP*EA; Vq = Vq1 + (EP * EA)), evaluate_if_numerical(V is V1+E). rev_eq_sum_product(P,A,F1,V,Vq):- rev(P,Q), rev(A,B), rev_eq_sum_product(Q,B,F,V,Vq), rev(F,F1). % % symbolic representation of sum. eqsum/2 and reqsum/2 % cited from: eba01.pl (28 May 2003) % ----------------------------------------------------------- % eqsum([],0). eqsum([X|Members],Sum):- eqsum(Members,Sum0), ( X=0 -> Sum = Sum0 ; ( Sum0=0 -> Sum = X ; Sum = Sum0 + X ) ). reqsum(A,B):- rev(A,C), eqsum(C,B). % % allocation % ----------------------------------------------------------- % allocation(N,A,[X|Y]):- allocation(N,A,A,[X|Y]). allocation(0,_,0,[]). allocation(N,A,B,[X|Y]):- integer(A), len([X|Y],N), allocation(_N1,A,B1,Y), % N1 is N - 1, % sum(Y,B1), K is A - B1 + 1, len(L,K), nth0(X,L,X), B is B1 + X. % % probability (percentile) by using allocation % ----------------------------------------------------------- % probabilities(0,[]). probabilities(N,[X|Y]):- integer(N), len([X|Y],N), allocation(N,100,[X|Y]). % % any ratio (weight) can be interpreted into a prob. % ----------------------------------------------------------- % scale(W,1/Z,P):- findall(Y,(kth(_K,W,X),Y is X/Z),P). probabilities(W,N,P):- len(W,N), sum(W,Z), scale(W,1/Z,P). % % degenerate probability %--------------------------------------------- degenerate(Ps):- kth(K,Ps,P), characteristic_vector(K,P,Ps,Ps). % % probability over base set with steps of levels. % ----------------------------------------------------------- % make_a_prob(P,base(M),steps(L)):- var(P), len(P,M), allocation(M,L,W), probabilities(W,M,P). make_a_prob(P,base(M),_):- \+ var(P), len(P,M), \+ ( member(P1,P), ( var(P1); P1 > 1; P1 < 0 ) ), sum(P,1). % % expected value % ----------------------------------------------------------- % expected_value(W,A,E/100):- len(A,N), probabilities(W,N,P), product_sum(P,A,_,E). % % expected value with equational % ----------------------------------------------------------- % expected_value_eq(W,A,E/100,Eq):- len(A,N), probabilities(W,N,P), product_sum_eq(P,A,_,E,Eq). % %--------------------------------------------------- % net present value (NPV) %--------------------------------------------------- % % time preference: discount factor %--------------------------------------------- interest_rate(1.1). discount_factor(R,Y,DF,DFV):- DF = R ^ (-Y), DFV is DF. npv(A,Y,Eq,V):- interest_rate(R), discount_factor(R,Y,DF,_), Eq = DF * A, V is Eq. % % conditional probabilities % ----------------------------------------------------------- % probability_of_event(W,E,P):- % conditionalization by event specified directly event(E), (var(E)->E = E1; sort(E,E1)), G = member(S,E1), findall(A,(probability(W,S,A),G),Ps), sum(Ps,P). probability_of_event(W,E,P,G):- \+ var(G), % conditionalization via constraints indirectly G=(Goal,M,[W,S,A]), % constraints with params findall([S1,A1], ( (M=do->(W=W1,S=S1,A=A1);true), probability(W1,S1,A1), Goal ), Xs), findall(S,member([S,A],Xs),E0), findall(A,member([S,A],Xs),Ps), sort(E0,E), sum(Ps,P). % % ----------------------------------------------------------- % % Utilities for list operations % ----------------------------------------------------------- % % edited: 14 Feb 2003. (imported from: set1.pl) % % index for tuples. % ----------------------------------------------------------- % % 1) only mention for a direct product of sets. index_of_tuple(B,A,Index):- \+ var(B), \+ var(A), len(B,LN), % base sets len(A,LN), len(Index,LN), findall(L, ( kth(K,B,BJ), %write(a(K,B,BJ)), kth(L,BJ,AJ),%write(b(L,BJ,AJ)), kth(K,A,AJ) %,write(c(K,A,AJ)),nl ), Index). index_of_tuple(B,A,Index):- \+ var(B), \+ var(Index), var(A), len(B,LN), % base sets len(Index,LN), len(A,LN), findall(AJ, ( kth(K,B,BJ), kth(K,Index,L), kth(L,BJ,AJ) ), A). % % descending/ascending natural number sequence less than N. % ----------------------------------------------------------- % dnum_seq([],N):-N<0,!. dnum_seq([0],1). dnum_seq([A|Q],N):- A is N - 1, len(Q,A), dnum_seq(Q,A). anum_seq(Aseq,N):-dnum_seq(Dseq,N),sort(Dseq,Aseq). dnum_seq1(Q,N):- M is N + 1, dnum_seq(Q0,M), subtract(Q0,[0],Q). anum_seq1(Q,N):- M is N + 1, anum_seq(Q0,M), subtract(Q0,[0],Q). % % inquire the goal multiplicity % ----------------------------------------------------------- % sea_multiple(Goal,Cond,N,M):- Clause=..Goal, findall(Cond,Clause,Z),len(Z,N),sort(Z,Q),len(Q,M). % bag0([],_A,0). bag0([C|B],A,N):- len([C|B],N), bag0(B,A,_N0), member(C,A). zeros(Zero,N):-bag0(Zero,[0],N). ones(One,N):-bag0(One,[1],N). % % subset_of/3 : subset-enumeration % ----------------------------------------------------------- % subset_of(A,N,As):- len(As,L), len(D,L), list_projection(D,As,B), len(B,N), sort(B,A). % a sequence of binary choice for a list: %-------------------------------------------------- list_projection([],[],[]). list_projection([X|Y],[_A|B],C):- X = 0, list_projection(Y,B,C). list_projection([X|Y],[A|B],[A|C]):- X = 1, list_projection(Y,B,C). % % characteristic_vector/3 % ----------------------------------------------------------- % % modified: 8 Feb 2003. without using nth1. % modified: 13 Feb 2003. without using member. characteristic_vector(X,B,Index):- \+ var(B), %member(X,B), list_projection(Index,B,[X]). characteristic_vector(1,X,[X|B],[1|DX]):- characteristic_vector(X,[X|B],[1|DX]). characteristic_vector(K,X,[_|B],[0|DX]):- characteristic_vector(K1,X,B,DX), K is K1 + 1. % % replace(Project,Goal,Base,Goal1):- % ----------------------------------------------------------- % % added: 15 Oct 2002. % a sequence of replacement of a subset of elements in Goal % which specified by a list, Project, 0-1^n, over Base % a list of length n, which brings about Goal1. % summary: % X=1 --> preserve the value of Base. % X=0 --> do replace with Goal1. replace([],[],[],[]). replace([X|A],[_|B],[Z|C],[Z|D]):- X = 0, replace(A,B,C,D). replace([X|A],[Y|B],[_|C],[Y|D]):- X = 1, replace(A,B,C,D). % % replace/4 % ----------------------------------------------------------- % % modified: 14 Feb 2003. bug fix. replace(K/N,L,S,L1):- \+ var(S), \+ var(L), len(L,N), len(L1,N), kth(K,L1,S), characteristic_vector(K,_S0,L,V), c_replace(V,L,L1,L1). % c_replace([],[],[],[]). c_replace([X|A],[_|B],[Z|C],[Z|D]):- X = 1, c_replace(A,B,C,D). c_replace([X|A],[Y|B],[_|C],[Y|D]):- X = 0, c_replace(A,B,C,D). % % sort without removal of duplicates %-------------------------------------------------- asort(A,B):- sort(A,C), bagof(CK, J^K^( kth(J,C,CK), kth(K,A,CK) ), B). % % permutation. % ----------------------------------------------------------- % % modified: 1 Sep 2002. to be used in is_neutral/2. % modified: 15 Oct 2002. add a non-variable constraint for the base set A. permutation_of(A,P,APs):- \+var(A), len(A,M), ordering(P,A,M), asc_nnseq(Qm,M), my_maplist(nth_of_permutation(A,P),Qm,APs). %caution: %notations 'A->B's bellow have cause precedence errors in IF-prolog. nth_poa(P,K,Ak->Pk):- alternatives(A), nth_of_permutation(A,P,K,Ak->Pk). nth_of_permutation(A,P,K,Ak->Pk):- len(A,M), ordering(P,A,M), nth_0(K,A,Ak), nth_0(K,P,Pk). permutate_of_order(PoA,[Q->R]):- poa(_P,PoA), alternatives(A), len(A,M), ordering(Q,A,M),%wn(Q), ordering(R,A,M),%wn(R), permutation(Q,PoA,R). permutation([],[],[]). permutation(Q,[A->P|PoA1],R):- subtract(Q,[A],Q1),nth_1(K,Q,A), subtract(R,[P],R1),nth_1(K,R,P), permutation(Q1,PoA1,R1). % % projection operator via index set. (an exchange economy?) % ----------------------------------------------------------- % % choice1 from base1 :: choice2 from base2. % modified : 15 Oct 2002. to be order-neutral (pending:-) % pcm([Choice1,Base1],[Choice2,Base2]):- pairwise_contract_map([Choice1,Base1],[Choice2,Base2]). pairwise_contract_map([Choice1,Base1],[Choice2,Base2]):- % len(Base1,N2), % len(Base2,N2), subset_of(Choice1,N1,Base1), subset_of(Choice2,N1,Base2), list_projection(Project,Base1,Choice1), list_projection(Project,Base2,Choice2). % % ppcm /2 using plist_projection % added: 15 Oct 2002. ppcm([Choice1,Base1],[Choice2,Base2]):- subset_of(C1,N1,Base1), subset_of(C2,N1,Base2), plist_projection(Project,Base1,Choice1,C1), plist_projection(Project,Base2,Choice2,C2). % % ----------------------------------------------------------- % % Utilities for outputs % ----------------------------------------------------------- % % redirect to file %-------------------------------- tell_goal(File,G):- (current_stream(File,write,S0)->close(S0);true), open(File,write,S), tell(File), nl, time_stamp('% file output start time ',_), nl, write('%---------- start from here ------------%'), nl, G, nl, write('%---------- end of data ------------%'), nl, time_stamp('% file output end time ',_), tell(user), close(S), % The following is to cope with the duplicated stream problem. (current_stream(File,write,S1)->close(S1);true). % 成功するゴールをすべて保存 %-------------------------------- tell_goal(File,forall,G):- G0 = (nl,write(G),write('.')), G1 = forall(G,G0), tell_goal(File,G1). tell_goal(File,forall_such_that,G,Condition):- % G should be occurred in the Condition. WRITE = (nl,write(G),write('.')), G1 = forall(Condition,WRITE), tell_goal(File,G1). % time stamp %-------------------------------- time_stamp(no,T):- get_time(U), convert_time(U,A,B,C,D,E,F,_G), T = [date(A/B/C), time(D:E:F)], nl. time_stamp(Word,T):- \+ var(Word), Word \= no, get_time(U), convert_time(U,A,B,C,D,E,F,_G), T = [date(A/B/C), time(D:E:F)], % format('~`.t~t~a~30|~`.t~t~70|~n',[Word]), write((Word,T)), nl. % len/2 : alternative to the length/2 len(A,B):- len_0(A,B). % length/2 of SWI-prolog fails at the case of both unbound variables. len_0(A,B):- var(A), var(B), !, fail. len_0([],0). len_0([_|A],N):- ((integer(N),N>0)->N0 is N -1; true), len_0(A,N0), ((integer(N0))->true; !,fail), N is N0 + 1. % kth/3 : alternative to the nth1/3 kth(K,Y,X):- kth_member(K,Y,X). kth_member(1,[X|Y],X):- \+ var(Y). kth_member(K,[_|Y],X):- \+ var(Y), kth_member(K1,Y,X), K is K1 + 1. % rev/2 : alternative to the reverse/2 rev(A,B):- len(A,L), len(B,L), rev(A,[],B), !. rev(A,A,[]):-!. rev(A,B,[C|D]):- rev(A,[C|B],D).