Why doesn't PostScript have arcus cosinus, and how to fix it?

Question

I found a formula for getting the angle between two 3D vectors (e.g. as shown in https://stackoverflow.com/a/12230083/6607497).

However when trying to implement it in PostScript I realized that PostScript lacks the arcus cosinus function needed to get the angle.

Why doesn't PostScript have that function and how would a work-around look like?

In Wikipedia I found the formula $\arccos(x)={\frac {\pi }{2}}-\arctan \left({\frac {x}{\sqrt {1-x^{2}}}}\right)}$ , but that looks a bit complicated; and if it's not: Why didn't they add acos (arcus cosinus) using that definition?

Did I miss something obvious?

Unless one of the original authors shows up, we are unlikely to get an answer regarding rationale. For resource efficiency, it was not uncommon for programming languages of the time (early 1980s) to restrict themselves to a minimal set of transcendental functions (`sin`, `cos`, {`log` | `ln`}, `exp`, `atan`) from which all others can be derived if needed. For a practical solution, you may want to look at [pst-math](https://www.ctan.org/ctan-ann/id/mailman.1348.1666106405.3715.ctan-ann@ctan.org) — njuffa, Oct 23 '22 at 21:19
@njuffa I had to dig, but in `pst-math.pro` I found: `/ACOS {dup dup mul neg 1 add sqrt exch atan DegToRad} bind def`. I probably don't want `DegToRad`. — U. Windl, Oct 23 '22 at 21:35
I am not convinced you need `acos`. In "How Futile are Mindless Assessments of Roundoff in Floating-Point Computation ?", W. Kahan suggests: ∠(x,y) := `2 * atan (norm (x*norm(y) - norm(x)*y) / norm (x*norm(y) + norm(x)*y)`. I have not tried it but it looks reasonable to use the Euclidean norm, and with Kahan we can be sure he knows what he is writing about. — njuffa, Oct 23 '22 at 21:50
@njuffa OK, but why do all the [math books](https://en.wikipedia.org/wiki/Euclidean_space#Angle) use `acos`? And can you formally define `norm(x)`? — U. Windl, Oct 23 '22 at 21:55
That is because math books deal with math, not some form of machine arithmetic. Many formulas found in math books are ill suited to numerical computation with limited precision. Kahan is the grand master of floating-point arithmetic, and has **deep** insights into how to compute *accurately* with it. In the cited document, Kahan actually uses the double-bar operator for the norm (but there is no MathJax available in this forum), and from context that is the Euclidean norm. It may be worthwhile to read Kahan's full discussion of angle computations in that paper. — njuffa, Oct 23 '22 at 22:01
I ran some experiments in C and Kahan's formula `2 * atan2 (norm (x * norm (y) - norm (x) * y), norm (x * norm (y) + norm (x) * y))` works well, as does the more commonly used `atan2 (norm (cross (x, y)), dot (x, y))`. Best I can tell, C's `atan2` is equivalent to Postscript's two-argument `atan` function. — njuffa, Oct 23 '22 at 23:21
See also [Numerically stable method for angle between 3D vectors](https://math.stackexchange.com/questions/1143354/numerically-stable-method-for-angle-between-3d-vectors) on Mathematics Stack Exchange. — njuffa, Oct 23 '22 at 23:39
This blog post about *accurate* angle computation also seems useful: https://www.jwwalker.com/pages/angle-between-vectors.html — njuffa, Oct 24 '22 at 00:48

score 2 · Accepted Answer · answered Oct 23 '22 at 21:50

As indicated in Wikipedia and done in pst-math, acos can be implemented rather easily using atan and sqrt (among other primitives) like this:

GS>% arcus cosinus, using degrees
GS>/acos { dup dup mul neg 1.0 add sqrt exch atan } bind def
GS>1 acos ==
0.0
GS>-1 acos ==
180.0
GS>0 acos ==
90.0

However this may be less efficient than a native implementation.

score 2 · Answer 2 · answered Oct 31 '22 at 10:08

As the comments indicated, if the only purpose for arcus cosinus is to get the angle between two vectors, another approach using the two-parameter atan is preferable.

I ran some experiments in C and Kahan's formula 2 * atan2 (norm (x * norm (y) - norm (x) * y), norm (x * norm (y) + norm (x) * y)) works well, as does the more commonly used atan2 (norm (cross (x, y)), dot (x, y)). Best I can tell, C's atan2 is equivalent to Postscript's two-argument atan function. – njuffa -Oct 23 at 23:21

So I tried to implement it (for 3D vectors (x y z)):

First the norm (length), named n3:

% norm 3D vector
%(x y z)
/n3 { dup mul exch dup mul add exch dup mul add sqrt } def

Then three helpers to multiply (v3m), add (v3a), and subtract (v3s) 3D vectors (for those not fearing stack rotations):

% multiply 3D vector with factor
%(x y z f)
/v3m {
    dup 3 1 roll mul 4 1 roll %(z*f x y f)
    dup 3 1 roll mul 3 1 roll %(z*f y*f x f)
    mul exch 3 -1 roll        %(x*f y*f z*f)
} def

% subtract 3D vectors
%(x1 y1 z1 x2 y2 z2)
/v3s {
    6 -3 roll               %(x2 y2 z2 x1 y1 z1)
    4 -1 roll sub           %(x2 y2 x1 y1 z1-z2)
    5 1 roll 3 -1 roll sub  %(z1-z2 x2 x1 y1-y2)
    4 1 roll exch sub       %(y1-y2 z1-z2 x1-x2)
    3 1 roll                %(x1-x2 y1-y2 z1-z2)
} def

% add 3D vectors
%(x1 y1 z1 x2 y2 z2)
/v3a {
    4 -1 roll add          %(x1 y1 x2 y2 z2+z1)
    5 1 roll 3 -1 roll add %(z2+z1 x1 x2 y2+y1)
    4 1 roll add           %(y2+y1 z2+z1 x1+x2)
    3 1 roll               %(x1+x2 y2+y1 z2+z1)
} def

Eventually the angle function (a) is a bit different, because I wanted to avoid heavy stack shuffling, so converted the vectors to arrays and assigned the arrays to a name each. I think it makes the code also a bit more readable, while a bit less efficient (you may try to use copy to duplicate the vectors' coordinates without using arrays and names (and the additional local dictionary):

% angle between two 3D vectors
%(x1 y1 z1 x2 y2 z2)
/a {
    [ 4 1 roll ]
    4 1 roll
    [ 4 1 roll ]
    %([v2] [v1])
    2 dict begin
    /v1 exch def
    /v2 exch def
        v1 aload pop n3 v2 aload pop
        4 -1 roll v3m
        %(|v1|*v2)
        v2 aload pop n3 v1 aload pop
        4 -1 roll v3m
        %(|v1|*v2 |v2|*v1)
        v3s n3
        %(||v1|*v2-|v2|*v1|)
        v1 aload pop n3 v2 aload pop
        4 -1 roll v3m
        %(|v1|*v2)
        v2 aload pop n3 v1 aload pop
        4 -1 roll v3m
        %(|v1|*v2 |v2|*v1)
        v3a n3
        %(||v1|*v2-|v2|*v1| ||v1|*v2+|v2|*v1|)
        atan
        %(atan(||v1|*v2-|v2|*v1|,||v1|*v2+|v2|*v1|))
        2.0 mul
    end
} def

Finally here are some test cases:

GS>0 0 1 0 0 1 a ==
0.0
GS>0 0 1 0 1 0 a ==
90.0
GS>0 0 1 0 1 1 a ==
45.0
GS>0 0 1 1 0 0 a ==
90.0
GS>0 0 1 1 0 1 a ==
45.0
GS>0 0 1 1 1 0 a ==
90.0
GS>0 0 1 1 1 1 a ==
54.7356071
%
GS>0 1 0 0 0 1 a ==
90.0
GS>0 1 0 0 1 0 a ==
0.0
GS>0 1 0 0 1 1 a ==
45.0
GS>0 1 0 1 0 0 a ==
90.0
GS>0 1 0 1 0 1 a ==
90.0
GS>0 1 0 1 1 0 a ==
45.0
GS>0 1 0 1 1 1 a ==
54.7356071
%
GS>0 1 1 0 0 1 a ==
45.0
GS>0 1 1 0 1 0 a ==
45.0
GS>0 1 1 0 1 1 a ==
0.0
GS>0 1 1 1 0 0 a ==
90.0
GS>0 1 1 1 0 1 a ==
59.9999962
GS>0 1 1 1 1 0 a ==
59.9999962
GS>0 1 1 1 1 1 a ==
35.264389
%
GS>1 0 0 0 0 1 a ==
90.0
GS>1 0 0 0 1 0 a ==
90.0
GS>1 0 0 0 1 1 a ==
90.0
GS>1 0 0 1 0 0 a ==
0.0
GS>1 0 0 1 0 1 a ==
45.0
GS>1 0 0 1 1 0 a ==
45.0
GS>1 0 0 1 1 1 a ==
54.7356071
%
GS>1 0 1 0 0 1 a ==
45.0
GS>1 0 1 1 1 0 a ==
59.9999962
GS>1 0 1 0 1 1 a ==
59.9999962
GS>1 0 1 1 0 0 a ==
45.0
GS>1 0 1 1 0 1 a ==
0.0
GS>1 0 1 1 1 0 a ==
59.9999962
GS>1 0 1 1 1 1 a ==
35.264389
%
GS>1 1 0 0 0 1 a ==
90.0
GS>1 1 0 0 1 0 a ==
45.0
GS>1 1 0 0 1 1 a ==
59.9999962
GS>1 1 0 1 0 0 a ==
45.0
GS>1 1 0 1 0 1 a ==
59.9999962
GS>1 1 0 1 1 0 a ==
0.0
GS>1 1 0 1 1 1 a ==
35.264389
%
GS>1 1 1 0 0 1 a ==
54.7356071
GS>1 1 1 0 1 0 a ==
54.7356071
GS>1 1 1 0 1 1 a ==
35.264389
GS>1 1 1 1 0 0 a ==
54.7356071
GS>1 1 1 1 0 1 a ==
35.264389
GS>1 1 1 1 1 0 a ==
35.264389
GS>1 1 1 1 1 1 a ==
0.0

I hope I didn't mess it up.

Nice job. There's a small function I use that helps make named local variables: `/args-begin { dup length dict begin {exch def} forall } def` To use it give it an array of variable names (in reverse stack order), eg. to use at the start of your `/a` function, it would look like: `/a { {z2 y2 x2 z1 y1 x1} args-begin ... `. Or to use in the middle where you're already using variables, it would look like: `{v1 v2} args-begin` and replaces 3 lines. — luser droog, Nov 03 '22 at 23:57

Why doesn't PostScript have arcus cosinus, and how to fix it?

2 Answers2