Our old Prolog system had a little computer algebra system (CAS) that supported vectors and marices. In the following we report on a new take that is less leaning towards object oriented programming and more leaning towards higher order logic programming.

Vectors will be represented as lists [A1,..,An]:

• Where A1 .. An are numbers.

If we would use index based access via a list predicate such as nth1/3, this would be slow. We demonstrate how to realize some standard canon of linear algebra operations with higher order list processing predicates only. We first look at the inner product:

    [W1,..,Wn] . [X1,..,Xn] = W1*X1 + .. + Wn*Xn

The above index based definition leads to a list based implementation where we first form the list of products [W1*X1,..,Wn*Xn] and then sum over the products. The Prolog code reads as follows:

:- ensure_loaded(library(lists)).

vecdot2(V, W, X) :-
    maplist(mul, V, W, H),
    foldl(add, H, 0, X).

mul(X, Y, Z) :- Z is X*Y.

add(X, Y, Z) :- Z is X+Y.

?- vecdot2([1,3,-5], [4,-2,-1], X).

We might want to avoid the materialization of the list H. It helps to view the sum in the form of a Horner schema, where we start with zero 0, and the subsequently form a product and a sum .. ((0+W1*X1)+W2*X2) .. +Wn*Xn. This leads to an alternative more direct implementation:

vecdot(V, W, X) :-
   foldl(muladd, V, W, 0, X).

muladd(X, Y, Z, T) :- T is X*Y+Z.

?- vecdot([1,3,-5], [4,-2,-1], X).

Matrixes will be represeted as lists [V1,..,Vm]:

• Where V1 .. Vm are vectors.

It is assume that the vectors V1 .. Vm have all the same size n, leading to a m x n matrix. The vectors are called the rows of the matrix, and the vector elements run along the column indexes. An operation that doesn't involve arithmetic is taking the transpose:

    (A^T)ij = Aji

Surprisingly we can use list processing to compute the transpose:

mattran([H|T], R) :- mattran(T, H, R).

mattran([], U, R) :- maplist(rdot([]), U, R).
mattran([J|T], H, R) :- mattran(T, J, TC), maplist(rdot, TC, H, R).

rdot(T, H, [H|T]).

?- mattran([[1,2],[3,4],[5,6]], X).

Linear algebra gets its name from linear equations of the form W1*X1 + .. + Wn*Xn = B. If they appear in plenty the weights Wij and targets Bi can be expressed as a matrix A respectively a vector y. And the equations can be expressed as a linear mapping A x = y.

matapply(M, V, W) :-
   maplist(vecdot(V), M, W).

?- matapply([[0,1],[1,0]], [2,3], X).

If there is a multiplicity of vectors x and y, we can arrange them into matrix B and C. The system of equations can then be expressed as matrix multiplication A B = C. Again list processing is enough to implement this operation:

matmul(M, N, R) :-
   mattran(N, H),
   maplist(matapply(H), M, R).

?- matmul([[1,2,3],[3,2,1]], [[4,5],[6,5],[4,6]], X).

We will turn to an application of linear algbra to number theory. It is possible to expess recurence relation such as the one of Fibonnaci numbers by a matrix. The initial condition and the recurrence relation is:

    F0 = 0,    F1 = 1,    Fn+2 = Fn+1 + Fn

We can use a vector Vn = [Fn+1, Fn]. The initial condition can now be expressed as V0 = [1, 0]. For the recurrence relation we need to find a matrix so that M * Vn = Vn+1. The following matrix does that:

    | 1 1 | | Fn+1 |   | Fn+1 + Fn |   | Fn+2 |
    |     | |      | = |           | = |      |
    | 1 0 | | Fn   |   |    Fn+1   |   | Fn+1 |

Lets check whether we can reproduce F10 = 55, F11 = 89 and F12 = 144:

?- matapply([[1,1],[1,0]], [89,55], X).

In general M^n V0 = Vn, a matrix exponentiation M^n can be implemented:

matexp(M, 1, M) :- !.
matexp(M, N, R) :- N mod 2 =:= 0, !,
   I is N // 2,
   matexp(M, I, H),
   matmul(H, H, R).
matexp(M, N, R) :-
   I is N-1,
   matexp(M, I, H),
   matmul(H, M, R).

We can again reproduce F12 = 144:

fib(N, X) :-
    matexp([[1,1],[1,0]], N, H),
    matapply(H, [1,0], [_,X]).

?- fib(12, X).

What are the last 8 digits of F1000000?

?- fib(1000000, _X), Y is _X mod 10^8.

The above will have a huge intermediate result _X, which the top-level doesn't show since we used a marked variable. This means we used bigint computation and only took modulo at the end. Alternatively we could have used modulo computation during the matrix exponentiation already.

Currently we only support numerical computation, either floats or bigint. We could implement Prolog operations that have strict vector or matrix signatures. Running it in the browser, we can answer in a blink "What are the last 8 digits of fib(1000000)?".