1:- module(matrix, [matrix_calc/3, matrix/2]).    2:- use_module(library(lists)).    3:- use_module(util(html)).    4:- use_module(util(misc)).    5:- use_module(util('meta2')).    6:- use_module(util('gb-lex-old')).    7:- use_module(pac(basic)).    8:- use_module(pac(reduce)).    9:- use_module(pac('expand-pac')).   10term_expansion --> pac:expand_pac.
   11
   12
   13:- use_module(pac(op)).   14:- use_module(pac(basic)).   15
   16equation_separator(;).
   17
   18% ?- qcompile(util(matrix)), module(matrix).
   19% ?- matrix:power_mod_table(3, X).
   20% ?- matrix:power_mod_table(3, X).
   21% ?- matrix:eqs_matrix((x+y=1;x-y=1), A, B, C).
   22% ?- matrix:eqs_matrix((x=1), A, B, C).
   23% ?- gblex:term_poly(x-1, _G3008).
   24
   25% ?- trace, matrix:matrix_calc([[[1, 2, 1], [1, 3, 1], [-2, -1, 0]], [[1, 2, 1], [1, 3, 1], [-2, -1, 0]], [[1, 0, 0], [0, 1, 0], [0, 0, 1]]], x+y, _3664).
   26
   27algebra(X, Y, Z) --> matrix(X, Y, Z).
   28%
   29is_matrix([_|_]).
   30
   31%
   32subst_rdiv(X, Y) :- subst_functor([((rdiv/2), (/))], X, Y).
   33
   34% ?- subst_functor([(f/1,  g)], f(g(f(1))), R).
   35% R = g(g(g(1))).
   36
   37subst_functor(_, X, X):- (number(X); var(X)), !.
   38subst_functor(S, X, Y):- X =..[F|A],
   39	maplist(subst_functor(S), A, B),
   40	functor(X, F, N),
   41	(  memberchk((F/N, G), S)
   42	->	H = G
   43	;	H = F
   44	),
   45	Y =..[H|B].
   46
   47% eval is necessary for @-expression, otherwise unnecessary.
   48% ?- disable_pac_query.
   49% ?- enable_pac_query.
   50% ?- matrix(inverse(square([1,2,3])),  R).
   51% ?- matrix(square([1,2,3]),  R).
   52% ?- eval(matrix:: (inverse @ square([1,2,3])),  R).
   53% ?- eval((:matrix):: (inverse @ square([1,2,3])),  R). % << Error.
   54% ?- eval(matrix@inverse@square@[1,2,3],  R).
   55% ?- eval(matrix@(inverse@(square@[1,2,3])),  R).
   56% ?- eval(matrix@(inverse@(square@[1,2,3])),  R), call(R, V).
   57% ?- eval(matrix:: unit(3), X).
   58% ?- eval(matrix:: (scalar(2+3)@unit(2)), X).
   59% ?- eval(matrix:: (unit(2) + scalar(2, unit(2))), X).
   60% ?- eval(matrix:: (unit(2) - scalar(2, unit(2))), X).
   61% ?- eval(matrix::(unit@ (:length([1,2]))), R).
   62% ?- eval(matrix::(unit(3) * unit(3)), X).
   63% ?- eval(matrix::(unit(3)+unit(3)), X).
   64% ?- eval(matrix::scalar(1+2, unit(3)+unit(3)), X).
   65% ?- eval(matrix::diag(square(vector(3, 1))), X).
   66% ?- eval(matrix::min(1,2), X).
   67% ?- eval(matrix::max(1,2), X).
   68% ?- matrix((2*unit(3))^1, X).
   69% ?- eval(matrix:: (2*unit(3))^1, X).
   70% ?- eval(matrix:: (2*unit(3))^2, X).
   71% ?- eval(matrix:: (2*unit(4))^2, X).
   72% ?- eval(matrix:: ((1+2)*unit(4))^2, X).
   73% ?- matrix((((1+2)*unit(4))^100), X), matrix(X * (X^(-1)), Y).
   74
   75:- bekind(matrix, []).   76square(L)	= :vector_to_square@L.
   77diag(M)		= :diagonal@M.
   78unit(N)		= :unit_matrix@N0 :- N0 is N.
   79inverse(M)	= xargs([M0]>>M0^(-1))@M.
   80vector(N, A)= V	:- length(V, N), maplist(=(A), V).
   81t(A)		= :transpose_matrix@A.
   82tr(A)		= :matrix_trace@A.
   83-A			= scalar(-1)@A.
   84A - B		= A + scalar(-1, B).
   85A + B		= :add_m_m@ A@B.
   86A * B		= :mul_m_m@ A@B.
   87s_add(A,B)	= :add_s_m(A0)@B :- A0 is A.
   88s_mull(A,B)	= :mul_s_m(A0)@B :- A0 is A.
   89A^N			= :exp_m(N0)@ A  :-  N0 is N.
   90max(A, B)	=  `C :-  ( A< B -> C=B; C=A).
   91min(A, B)	=  `C :-  ( A< B -> C=A; C=B).
   92scalar(A, X)= :mul_s_m(A0)@ X :- A0 is A.
   93subst(X, P)	= :poly_subst@X@ (math:poly(P)).
   94X			= `X.
   95:- ekind.
   96
   97% ?- matrix_calc([[[1,2,1],[1,3,1],[-2,-1,0]],[[1,2,1],[1,3,1],[-2,-1,0]],[[1,0,0],[0,1,0],[0,0,1]]], x+y, R).
   98% ?- matrix_calc([[[1,2,1],[1,3,1],[-2,-1,0]],[[1,2,1],[1,3,1],[-2,-1,0]],[[1,0,0],[0,1,0],[0,0,1]]], x(1)+x(2), R).
   99
  100:- bekind(matrix_calc//1,[]).
  101x		= with([A|_], `A).
  102y		= with([_,A|_], `A).
  103z		= with([_,_,A|_], `A).
  104x(K)		= with(L, :nth1(K, L)).
  105square(L)	= :vector_to_square@L.
  106diag(M)		= :diagonal@M.
  107unit(N)		= :unit_matrix@N0 :- N0 is N.
  108inverse(M)	= xargs([M0] >> M0^(-1))@M.
  109vector(N, A)= V	:- length(V, N), maplist(=(A), V).
  110t(A)		= :transpose_matrix@A.
  111tr(A)		= :matrix_trace@A.
  112-A			= scalar(-1)@A.
  113A - B		= A + scalar(-1, B).
  114A + B		= :add_m_m@ A@B.
  115A * B		= :mul_m_m@ A@B.
  116s_add(A,B)	= :add_s_m(A0)@B :- A0 is A.
  117s_mull(A,B)	= :mul_s_m(A0)@B :- A0 is A.
  118A^N			= :exp_m(N0)@ A  :-  N0 is N.
  119max(A, B)	= `C :-  ( A< B -> C=B; C=A).
  120min(A, B)	= `C :-  ( A< B -> C=A; C=B).
  121scalar(A, X)= :mul_s_m(A0)@ X :- A0 is A.
  122subst(X, P)	= :poly_subst@ X @ (math:poly(P)).
  123X