# Manipulating matrices for 3D drawing

I've written a short library for manipulating matrices for 3D drawing. I've tried to strike a balance between speed and readability. Anything to improve?

%!
%mat.ps
%Matrix and Vector math routines

/ident { 1 dict begin /n exch def
[
1 1 n { % [ i
[ exch % [ [ i
1 1 n { % [ [ i j
1 index eq { 1 }{ 0 } ifelse % [ [ i b
exch % [ [ b i
} for % [ [ b+ i
pop ] % [ [ b+ ]
} for % [ [b+]+ ]
]
end } def

/rotx { 1 dict begin /t exch def
[ [ 1  0      0         ]
[ 0  t cos  t sin neg ]
[ 0  t sin  t cos     ] ]
end } def

/roty { 1 dict begin /t exch def
[ [ t cos      0  t sin ]
[ 0          1  0     ]
[ t sin neg  0  t cos ] ]
end } def

/rotz { 1 dict begin /t exch def
[ [ t cos  t sin neg  0 ]
[ t sin  t cos      0 ]
[ 0      0          1 ] ]
end } def

/.error where { pop /signalerror { .error } def } if

/dot { % u v
2 copy length exch length ne {
/dot cvx /undefinedresult signalerror } if
% u v
0 % u v sum
0 1 3 index length 1 sub { % u v sum i
3 index 1 index get exch % u v sum u_i i
3 index exch get % u v sum u_i v_i
mul add % u v sum
} for % u v sum
3 1 roll pop pop % sum
} bind def

% [ x1 x2 x3 ] [ y1 y2 y3 ]  cross  [ x2*y3-y2*x3 x3*y1-x1*y3 x1*y2-x2*y1 ]
/cross { % u v
dup length 3 ne 2 index length 3 ne or {
/cross cvx /undefinedresult signalerror } if
% u v
exch aload pop 4 3 roll aload pop % x1 x2 x3 y1 y2 y3
[
5 index 2 index mul % ... [ x2*y3
3 index 6 index mul sub % ... [ x2*y3-y2*x3
5 index 5 index mul % ... [ x2*y3-y2*x3 x3*y1
8 index 4 index mul sub % ... [ x2*y3-y2*x3 x3*y1-x1*y3
8 index 5 index mul % ... [ x2*y3-y2*x3 x3*y1-x1*y3 x1*y2
8 index 7 index mul sub % ... [ x2*y3-y2*x3 x3*y1-x1*y3 x1*y2-x2*y1
]
7 1 roll 6 { pop } repeat
} bind def

/transpose { STATICDICT begin
/A exch def
/M A length def
/N A 0 get length def
[
0 1 N 1 sub { /n exch def
[
0 1 M 1 sub { /m exch def
A m get n get
} for
]
} for
]
end } dup 0 6 dict put def

/matmul { STATICDICT begin
/B exch def
B 0 get type /arraytype ne { /B [B] def } if
/A exch def
A 0 get type /arraytype ne { /A [A] def } if
/Q B length def
/R B 0 get length def
/P A length def
Q A 0 get length ne {
/A A transpose def
/P A length def
Q A 0 get length ne {
A B end /matmul cvx /undefinedresult signalerror
} if
} if

[
0 1 P 1 sub { /p exch def % rows of A
[
0 1 R 1 sub { /r exch def % cols of B
0
0 1 Q 1 sub { /q exch def % terms of sum
A p get q get
B q get r get mul
} for
} for
]
} for
]

end } dup 0 10 dict put def

%u v {operator}  vop  u(op)v
%apply a binary operator to corresponding elements
%in two vectors producing a third vector as result
/vop { 1 dict begin
/op exch def
2 copy length exch length ne {
/vop cvx end /undefinedresult signalerror
} if

[ 3 1 roll % [ u v
0 1 2 index length 1 sub { % [ ... u v i
3 copy exch pop get % u v i u_i
3 copy pop get      % u v i u_i v_i
op exch pop         % u v u_i(op)v_i
3 1 roll            % u_i(op)v_i u v
} for % [ ... u v
pop pop ]

end } def

%length of a vector
/mag { 0 exch { dup mul add } forall } def

/elem { % M i j
3 1 roll get exch get % M_i_j
} def

/det {
dup length 1 index 0 get length ne { /det cvx /typecheck signalerror } if
1 dict begin /M exch def
M length 2 eq {
M 0 0 elem
M 1 1 elem mul
M 0 1 elem
M 1 0 elem mul sub
}{
M length 3 eq {
}{ /det cvx /rangecheck signalerror } ifelse
} ifelse
end
} def

-
Current version: github.com/luser-dr00g/xpost/blob/master/data/mat.ps – luser droog Sep 16 '15 at 21:31
This code is used to produce the images in these posts across the SE network: 1 2 3 – luser droog Jan 30 at 7:25
– luser droog Jan 30 at 7:25

As far as readability goes, maybe you'd like to add an arg1 arg2 arg3 -- result1 result2 type of comment to each function? So we know what was on the stack before invoking, and what we expect to find afterward? Almost equivalently, you might show a unit test invoking your function and verifying the result.

-
Good advice. I usually do that and wonder why I didn't here. I'm working up a new interpreter to run it, but it will have these things before entering the repository. :) – luser droog May 2 '13 at 6:42
Oh, I remember! I was hoping the stylized local argument convention would make the parameters apparent on the functions without the usual comment. They either construct the argument directly where the layout naturally illustrates the result, or contain intermediate stack-comments up through the final line, illustrating the result. But, it appears it's too clever which is bad. The entire lack of comments in the determinant function det was an oversight. Rats! matmul and transpose need comments, too. :( – luser droog May 2 '13 at 6:49

I just noticed some issues with the vop procedure. It's not nestable the way it's written. So you can't use it to add two matrices like m1 m2 { {add}vop }vop. I've got a favorite trick using token that'll fix this. But I've got serious versioning issues, because I've post this code everywhere with no identifying marks!! I'm an idiot.

So another big improvement would be, maybe, a copyright notice, um, a release date, a version number?!

I've also got a handful of routines that I've considered adding because I keep needing them in the same programs that use this. But they're not really related to matrices, so I've kept them out to avoid the kitchen-sink problem.

/div { dup 0 eq { pop 10000 }{ div } ifelse } bind def
/normalize{[1 index mag 2 index length 1 sub{dup}repeat]{div}vop}def
/tan { dup sin exch cos div } def


Replacing div with a version that doesn't explode on divide-by-zero has gotten things to run that otherwise wouldn't, often enabling me to locate the problem visually.

normalize probably should be in there, if it weren't for the versioning issue.

And tan is curiously omitted from the PostScript operator set. Not sure if this is the best way, but it's probably the simplest. But this requires the div guard, because tangents go to infinity when the cos is 0. The practical selection of infinity is perhaps best left to the individual program. So, I suppose those two should stay out.

-

Much later, I notice an error in the matmul routine if the right argument is a vector, it should be promoted to a column vector in order to multiply with a left matrix. This code doesn't do that.

Rather than

B 0 get type /arraytype ne { /B [B] def } if


it needs to make unit subarrays of each element of B.

B 0 get type /arraytype ne {
/B [ 0 1 B length 1 sub {
B exch 1 getinterval
} for ] def
} if


Oh, and mag should probably take a sqrt.

And det should signal the undefinedresult error rather than typecheck if the argument is not a square matrix. Perhaps it should signal a typecheck is the argument is not an array or scalar number.

-