Matrices

[benryves]

I recently received an email from Tom Lake requesting the addition of matrix functions to BBC BASIC. I must confess that I'd never seen a BASIC dialect that offered native matrix functions, but some digging located the following:

• PowerBASIC MAT
• HP BASIC MAT =/INPUT/PRINT/READ
• Wang BASIC MAT

I'm wondering if people had any requirements or suggestions for the MAT statement. There are some limitations I have to work with; for example, any matrix operation must be prefixed with MAT, and matrices are not efficiently resizeable as BBC BASIC only stores the current size, not the maximum size (matrices could be shrunk, but not expanded, without leaking memory).

To whet your appetites, here's the start; I've written some parsing code for the simplest MAT statements - ZER (set all elements to zero), CON (set all elements to one) and CON(expr) (set all elements to expr). Though I've just realised that I'm missing the = between the result array and the expression to the right, whoops.

[tr1p1ea]

Wow, matrix manipulation? I am very impressed!

[benryves]

As well as mathematical uses I'm sure that the ability to quickly copy one array's contents to another or reset all the elements to a particular value would be quite useful.

I'd like to try and be fairly practical with this, but am not sure of the most sensible way to do it. Would multiplying a two dimensional array ("matrix") with a one-dimensional array ("vector") transform the vector by the matrix? Would a()^b() calculate the cross product and a()*b() the dot product if both were single-dimension arrays, or would you use a()*b() and a().b() respectively (or maybe even prefix those statements with VEC instead of MAT)?

I'm concerned about doing this sensibly, and not ending up releasing something with a silly/non-standard syntax that I need to change at a later date.

[toml_12953]

I'm wondering if people had any requirements or suggestions for the MAT statement. There are some limitations I have to work with; for example, any matrix operation must be prefixed with MAT, and matrices are not efficiently resizeable as BBC BASIC only stores the current size, not the maximum size (matrices could be shrunk, but not expanded, without leaking memory).

Why not just go with one of the implementations you list? Those are very good (especially the Wang) One thing, though. the operation of the determinant function, DET, is slightly different in different compilers. In some, DET takes no parameters. It only retuns the determinant of the most recently inverted matrix. In order to find the determinant of array A, you'd need to do this

Code: Select all

MAT T = INV(A)
PRINT DET

In others, DET takes the name of an array and there's no need to invert the matrix first. You can do this:

Code: Select all

PRINT DET(A)

Some compilers allow both forms. The ANSI/ISO standard requires DET to have an argument. I like having a choice but if you only have room for one, I'd go with the ISO standard, DET(A).

Tom Lake

[benryves]

Thanks for the suggestion. Unfortunately, I would not be able to add a DET function in that form - I can only add statements to BBC BASIC, not functions. To that end, the former method would work, as I could automatically create/modify a DET variable after inverting a matrix - but this would slow things down unnecessarily if you wished to invert a matrix but didn't want to know the determinant.

The following could work:

Code: Select all

MAT d=DET(a())

however, you would not be able to use DET() as a regular function, eg in a PRINT statement.

[toml_12953]

The following could work:

Code: Select all

MAT d=DET(a())

however, you would not be able to use DET() as a regular function, eg in a PRINT statement.

But DET is a single number, not a matrix. It's used to determine how "good" an inversion is. The larger DET is, the better the inversion; the less loss of precision. If DET is zero, then either the matrix is singular (has no inverse) or its inverse is useless due to lost precision.

Code: Select all

MAT d=DET(a())

isn't conceptually right at all since "MAT d" should always refer to a matrix named d, not a scalar. Since you have to calculate the determinant during the inversion of a matrix anyway, there shouldn't be too much overhead to assign that value to the system variable DET (with no parameter). It's done this way in Wang BASIC and many others.

Tom Lake

[benryves]

Point taken. I shall use the parameterless DET method, then.

What is your view on the different dimensions of arrays? For some operations it wouldn't really matter, I suppose (ZER or CON) and with some it would be nonsensical (assigning a one-dimensional array to a two-dimensional array), but what would you expect to happen if, say, you multiplied a one-dimensional array with a two-dimensional one?

[toml_12953]

Point taken. I shall use the parameterless DET method, then.

What is your view on the different dimensions of arrays? For some operations it wouldn't really matter, I suppose (ZER or CON) and with some it would be nonsensical (assigning a one-dimensional array to a two-dimensional array), but what would you expect to happen if, say, you multiplied a one-dimensional array with a two-dimensional one?

The standards and conventions are very clear on this. Here's a quote from a BASIC manual:

If one array operand of the multiplication operator is one-dimensional and the other is two-dimensional, the product will be one-dimensional. If the first array operand is one-dimensional, it is treated as a row vector (single row with multiple columns) and must match the first dimension of the second array. If the second array operand is onedimensional, it is treated as a column vector (many rows in one column) and must match the second dimension of the first array.
In other words an array with m elements may be multiplied by an m by n array, or an L by m array may be multiplied by an array with m elements. In the first case the product will be a one-dimensional array with n elements, while in the second case the product will be a one-dimensional array with L elements.

[benryves]

You can download a preliminary version of the matrix additions by following this link.

Currently implemented:

• MAT b()=a() (assignment/copying)
• MAT a()=ZER (set all elements to zero/empty string)
• MAT a()=CON (set all elements to one)
• MAT a()=CON(expr) (set all elements to result of expr)
• MAT a()=IDN (set to identity matrix)
• MAT c()=a()+b() (add all elements)
• MAT c()=a()-b() (subtract all elements)
• MAT c()=a()*b()
• MAT b()=TRN(b()) (transposition)
• MAT PRINT a() (display a matrix's contents)

To do:

• MAT b()=INV(a()) (inversion)
• DET (determinant, set by INV())
• MAT c=a()*(scalar)
• MAT READ (read array elements from DATA statements)
• MAT INPUT (input array elements from prompt)

I'm not sure how many features I'll be able to fit in (I'm running out of ROM space) but I'll do my best.

I'm sure there will be bugs in the features I've implemented above, so please let me know when you find them!

Edit: Changed release 749 to 750; this includes a quick hack to allow multiplication of single-dimensional matrices, which hasn't been documented but should work.

[toml_12953]

I'm sure there will be bugs in the features I've implemented above, so please let me know when you find them!

I don't know if you're planning on implementing this or not but

Code: Select all

MAT PRINT A;

should print matrix A in a packed format while

Code: Select all

MAT PRINT A (or A,)

should print them in the usual format (just like a regular PRINT statement)

Also, are the parens really necessary?

Code: Select all

MAT A=B

should be sufficient to the interpreter.

Tom L

[benryves]

The parentheses are optional wherever an array is specified, but I've included them in all of the examples for clarity (the documentation is still incomplete).

Here is build 760. I've added a few bug fixes, including an accuracy fix on matrix multiplication (it used to store the intermediate values in the result, which resulted in very poor results when the result matrix was an integer array) and a crashing bug fix when array names are zero characters long.

I have also added the compacted form of MAT PRINT a; - is there any accepted way to suppress the line feed between PRINTed matrices?

I have also modified quite a lot of the code to try and gain a little more space; I'm down to 490 bytes to squeeze the inversion/determinant code in, which is more than a little tight!

[toml_12953]

I have also added the compacted form of MAT PRINT a; - is there any accepted way to suppress the line feed between PRINTed matrices?

Not in standard BASIC. Maybe you can eliminate the check for parens when MAT is used. That might free up some space.

When I tried this:

Code: Select all

10 DIM A(2,2)
20 MAT A=IDN
30 MAT PRINT A;

I get

011
111
111

instead of the correct

100
010
001


Tom Lake

[benryves]

Maybe you can eliminate the check for parens when MAT is used. That might free up some space.

It would free some space, but not a considerable amount. I'll see what I can do.

When I tried this:

Code: Select all

10 DIM A(2,2)
20 MAT A=IDN
30 MAT PRINT A;

I get

011
111
111

instead of the correct

100
010
001

Whoops, sorry. Build 761 fixes that (a rather overzealous optimisation).

[toml_12953]

I have also modified quite a lot of the code to try and gain a little more space; I'm down to 490 bytes to squeeze the inversion/determinant code in, which is more than a little tight!

Maybe you could eliminate Revision 752
"Added check to ensure that result matrix is neither of the operand matrices when multiplying."

This should be legal (and is in standard BASIC).

Code: Select all

MAT X = A * X

Should be allowed since the right side should be calculated then the result assigned to the matrix on the left of the equal sign (Unless you're multiplying in place to save space!)

Tom Lake

[benryves]

Yes, matrix operations are currently done in-place. This is the same reason that you cannot transpose a matrix onto itself. I have very limited access to the internals of BBC BASIC as far as variable allocation, variable manipulation and parsing are related, which means I have to implement these features myself. As (for example) I haven't written a string allocation routine yet either you can't copy an array of strings at all. This is also why all arrays have to be dimensioned beforehand.

I'd like to make the implementation as sensible as possible, don't worry, but these things take time (something I don't have a lot of at the moment!)

I don't suppose you have any recommendations or suggestions for an algorithm that could be used to calculate the inversion of a matrix?

Thanks for all of your help and suggestions!