Take the 2-minute tour ×
Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.

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
                add
            } 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 {
                M aload pop cross dot
            }{ /det cvx /rangecheck signalerror } ifelse
        } ifelse
    end
} def
share|improve this question
add comment

2 Answers

up vote 2 down vote accepted

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.

share|improve this answer
    
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
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.