This shows that numerical analysis algorithms can be implemented elegantly in Prolog.
check_determ domains scalar = real vector = scalar* matrix = vector* predicates gauss_jordan_elimination(matrix,matrix) gauss_jordan_elimination_1(matrix,matrix,matrix) gaussian_elimination_with_backsubstitution(matrix,vector) gaussian_elimination(matrix,matrix,matrix) backsubstitution(matrix,vector,vector) backsubstitution_1(vector,vector,scalar,scalar) pivot_row(matrix,vector,matrix) normalized_pivot_row(vector,vector) normalized_pivot_row_1(vector,scalar,vector) normalized_other_rows(matrix,vector,matrix) normalized_other_row(vector,scalar,vector,vector) reverse(matrix,matrix,matrix) clauses /*-------------------------------------------------------------------------*/ /* Gauss-Jordan Elimination with Partial Pivoting */ /* It handles any number of right-hand sides */ /* It can be used to invert a matrix */ /* */ /* e.g. gauss_jordan_elimination([[2,1,3, 13, 7], Gives [[1,1], */ /* [3,2,1, 10, 8], [2,2], */ /* [2,4,1, 13,11]], X). [3,1]] */ /*-------------------------------------------------------------------------*/ gauss_jordan_elimination(Ass, Xss):- gauss_jordan_elimination_1(Ass, [], Xss). gauss_jordan_elimination_1([], ReversedXss, Xss):- !, reverse(ReversedXss, [], Xss). gauss_jordan_elimination_1(Lower, Upper, Xss):- pivot_row(Lower, PivotRow, OtherRows), normalized_pivot_row(PivotRow, NormalizedPivotRow), normalized_other_rows(OtherRows, NormalizedPivotRow, NewLower), normalized_other_rows(Upper, NormalizedPivotRow, NewUpper), gauss_jordan_elimination_1(NewLower, [NormalizedPivotRow|NewUpper], Xss). pivot_row([PivotRow], PivotRow, []). pivot_row([[X1|Row1],[X2|Row2]|Rows], PivotRow, [[X2|Row2]|Rest]):- abs(X1) > abs(X2), !, pivot_row([[X1|Row1]|Rows], PivotRow, Rest). pivot_row([Row1,Row2|Rows], PivotRow, [Row1|Rest]):- pivot_row([Row2|Rows], PivotRow, Rest). normalized_pivot_row([0|_], _):- !, write("No non-zero pivot\n"), fail. normalized_pivot_row([Pivot|As], Bs):- normalized_pivot_row_1(As, Pivot, Bs). normalized_pivot_row_1([], _, []). normalized_pivot_row_1([A|As], Pivot, [B|Bs]):- B = A / Pivot, normalized_pivot_row_1(As, Pivot, Bs). normalized_other_rows([], _, []). normalized_other_rows([[A|As]|Ass], PivotRow, [Bs|Bss]):- normalized_other_row(As, A, PivotRow, Bs), normalized_other_rows(Ass, PivotRow, Bss). normalized_other_row([], _, [], []). normalized_other_row([A|As], X, [P|Ps], [B|Bs]):- B = A - X * P, normalized_other_row(As, X, Ps, Bs). /* reverse(Xs, Ys, Zs) is true if Zs is the result of concatenating the */ /* reverse of the list Xs to the front of the list Ys. */ reverse([], Ys, Ys). reverse([X|Xs], Ys, Zs):-reverse(Xs, [X|Ys], Zs). /*-------------------------------------------------------------------------*/ /* Gaussian Elimination with Partial Pivoting and Backsubstitution */ /* It handles only one right-hand side */ /* */ /* e.g. gaussian_elimination_with_backsubstitution([[2,1,3, 13], gives [1, */ /* [3,2,1, 10], 2, */ /* [2,4,1, 13]], X). 3] */ /*-------------------------------------------------------------------------*/ gaussian_elimination_with_backsubstitution(Ass, Xs):- gaussian_elimination(Ass, [], ReversedReducedAss), backsubstitution(ReversedReducedAss, [], Xs). gaussian_elimination([], Xss, Xss):-!. gaussian_elimination(Lower, Upper, Xss):- pivot_row(Lower, PivotRow, OtherRows), normalized_pivot_row(PivotRow, NormalizedPivotRow), normalized_other_rows(OtherRows, NormalizedPivotRow, NewLower), gaussian_elimination(NewLower, [NormalizedPivotRow|Upper], Xss). backsubstitution([], Xs, Xs). backsubstitution([As|Ass], Ys, Xs):- backsubstitution_1(Ys, As, 0, Y), backsubstitution(Ass, [Y|Ys], Xs). backsubstitution_1([], [B], Acc0, Acc):- Acc = B - Acc0. backsubstitution_1([X|Xs], [A|As], Acc0, Acc):- Acc1 = Acc0 + X * A, backsubstitution_1(Xs, As, Acc1, Acc).