[cairo] Shouldn't Cairo use/offer degrees rather than radians?

[David Kastrup]

Bill Spitzak <spitzak at gmail.com> writes:

> Your matrix for 45 degrees has equal values in all 4 entries. When
> this matrix is mulitplied by itself you get a matrix made of
> 0.9999999999999998 and 0.0. (That value is actually
> 1.0-2.220446049250313e-16).
>
> If you use the normal sin and cos values (which are unequal to each
> other)

Because M_PI/4 is not an exact representation of pi/4 (when using the
standard i386 FPU, you'd have to use 80 bits of precision rather than
64 bits anyway).  I seem to remember that the FPU actually has a command
for generating its own version of pi, FLDPI.  But the compiler does not
convert uses of M_PI into it.

> you get 1.0 and 2.220446049250313e-16 which is the same magnitude of
> error.

The same total _absolute_ error per element.  Getting 2e-16 instead of 0
is a larger relative error than getting (1-2e-16) instead of 1.

> However squaring 0.9999999999999998 gets a number further away from
> 1.0, while squaring 2.220446049250313e-16 goes closer to zero. For
> many operations (such as further multiplies, getting the determinant,
> etc) this means the second matrix is more accurate.

The determinant is a _square_ measure.  The actual scale factor for
geometric operations is its root.  Further multiplies will increase the
inaccuracy.  Here is how this works: we start with 45 degrees being off
as

(sqrt(0.5)+eps -sqrt(0.5)+eps
 sqrt(0.5)-eps sqrt(0.5)+eps)

After squaring (omitting larger powers of eps itself) we have

(sqrt(8)*eps -1
 1   sqrt(8)*eps)

Square again and we get

(-1  -sqrt(32)*eps
 sqrt(32)*eps  - 1)

Square a final time (we would want to get unity again) and we arrive at

(1 sqrt(128)*eps
 -sqrt(128)*eps 1)

Each further squaring will double the error terms (as long as we can
ignore eps^2 terms).

Ok.  Now let's start with the situation of my patch which can be
described as

(sqrt(0.5)+eps -sqrt(0.5)-eps
 sqrt(0.5)+eps sqrt(0.5)+eps)

Squaring once gives us

(0 -1-sqrt(8)*eps
 1+sqrt(8)*eps 0)

Squaring twice another time gives

(-1-sqrt(32)*eps   0
 0 -1-sqrt(32)*eps 0)

Another time gives

(1+sqrt(128)*eps 0
 0 1+sqrt(128)*eps)

Each further squaring will double the error terms (as long as we can
ignore eps^2 terms).  Since the magnitude of the term in question is
dominated by 1 (rather than by some factor of eps), eps^2 terms have
less opportunity to accumulate than in the sin/cos case.  So we get the
same kind of linear error growth as with sin/cos, a somewhat better
start of quadratic error growth, and perfect angles.

The initial eps in the latter case may start out larger, but your data
does not even show that.  The subsequent growth will be the same until
terms of eps^2 become relevant which is sooner the case for sin/cos than
for my patch.  The angle remains perfect for the matrix generated by my
patch while its error grows for the sin/cos matrix in use currently.

>>> You wanted two rotates by 45 degrees to be perfect.
>>
>> Uh, no?  Two rotates by 45 degrees illustrate a _compromise_.  The
>> degree of rotation will generally be a perfect 90 degrees because of all
>> rotation matrix elements having the same magnitude, the total magnitude
>> (determinant of the scaling matrix I think) will lightly be slightly
>> more wrong than the magnitude of the radian API (which likely also fails
>> to be 1 due to sin/cos of M_PI/4 in floating point also being
>> approximations).
>
> As I tried to show above the answer current Cairo gets for pi/4 is
> better than the one you get for 45 degrees.

If by "better" you mean that you ignore one side of the compromise,
namely that executing two 45 degree turns in sequence results in an
_exact_ angle of 90 degrees.

> Your function could be improved by using sin and cos rather than
> trying to do sin(pi/2-a) in place of cos in order to make sin and cos
> equal for 45 degrees.

I fail to see a net advantage.

>> Did you even read the patch and its rationale?  Or are you making up
>> that straw man on the fly?  Multiples of 90 degrees are perfect.
>> There are currently several fast paths in Cairo's code paths which
>> actually _check_ for that kind of perfection.
>
> I agree there are fast paths that are not allowing such errors in the
> matrix, and that changing cairo_rotate to detect angles near n*pi/2
> and produce 1.0 and 0.0 would probably be a much faster fix than
> trying to track down all the fast path mistakes. It would also remove
> suprises in the matrix for users.

"detect angles near n*pi/2" will introduce discontinuities.  My patch
does not introduce discontinuities because it changes functions in
longer stretches of 1.0.

>>> The exact same value is stored whether the rotation is sent as degrees
>>> or radians.
>>
>> DID YOU LOOK AT THE PATCH?????  I cannot believe you did when you state
>> this.
>
> Yes your patch in effect finds the quadrant the angle is in. It
> substitutes sin(pi/2-x) for cos(x) in order to make a right angle have
> a cos of 0.0, however as stated above I believe this is a mistake as
> the determinant of the matrix is not 1.0 for other angles.

It isn't 1.0 for other angles in general currently either.

> My version instead finds the closest axis (or the quadrant if you
> rotate by 45 degrees). The angle passed to sin/cos is in the range
> -pi/4..pi/4. This allows the normal sin/cos functions to be used.

Your quadrant switch is at a point where it will cause a detectable
discontinuity.

> Here it is (in python sorry):
>
> from math import *
> M_PI_2 = pi/2

You are working in radians: this does not work since multiples of pi/2
do not have a reliable representation in floating point arithmetic.
I may be repeating myself here.  Putting your unwillingness to
acknowledge the principal point of this thread aside, the basic idea of
your code can of course be implemented in the same manner in degrees.

It was my first approach as well (and actually, I admit it _was_ done in
radians at first too).  The main problem with radians as opposed to
degrees is that "fromaxis" is not reliably zero for multiples of 90
degrees since multiples of M_PI/2 have no dependable numeric
representation.

> # this is needed to avoid producing negative zero:
> def neg(x):
>   return -x if x else 0.0
>
> def sincos(a):
>   axis = round(a / M_PI_2)
>   fromaxis = a - axis * M_PI_2
>   s = sin(fromaxis)
>   c = cos(fromaxis)
>   x = int(axis) & 3
>   if   (x==0): return(     s,     c)
>   elif (x==1): return(     c, neg(s))
>   elif (x==2): return(neg(s), neg(c))
>   else:        return(neg(c), s)
>
> If you think degrees would help then substitute 90 for M_PI_2 and put
> fromaxis*(M_PI_2/90) to the sin and cos function calls. Though
> intuitively it seems like 90*round(a/90) is somehow more accurate than
> M_PI_2*round(a/M_PI_2) I have not been able to find a value where this
> actually happens so I don't believe the degree version is necessary.

Belief is one thing, verification another.  For example, I have

#include <math.h>
#include <stdio.h>

int
main()
{
  double da = 990.0;
  double ra = 11*(M_PI_2);
  printf ("%g %g\n", da-90.0*round(da/90.0),
	  ra-M_PI_2*round(ra/M_PI_2));
  return 0;
}

And get as result

-*- mode: compilation; default-directory: "/tmp/" -*-
Compilation started at Wed Jul 12 10:44:18

gcc gaga.c -lm && ./a.out
0 -1.77636e-15

Compilation finished at Wed Jul 12 10:44:18

Now I'll readily admit that 11 was more than I expected here: the
problem is apparently masked through the rules of IEEE arithmetic more
often than I would have expected.

>>> 2. Add an api to rotate so the x axis passes through an x,y point (and
>>> possibly scales by hypot(x,y)). This would provide an "accurate"
>>> rotate for a whole lot of cases that actually come up in real
>>> graphics.
>>
>> For multiples of 90 degrees, you can "trivially" just specify the
>> transform matrix, yet nobody does.  This is not how people think, and we
>> are talking about an API for people.
>
> I propose cairo_rotate_atan(x)

I don't think a 2-quadrant version is a good idea at all.  If you want
that, you can use cairo_rotate_atan2(x,1) explicitly,

> and cairo_rotate_atan2(x,y) which are exactly the same as
> cairo_rotate(atan(x)) and cairo_rotate(atan2(x,y)) except for greater
> accuracy.

By the way, "x,y" instead "y,x" for atan2 is _really_ asking for
confusion.

"Greater accuracy" again is a euphemism: we are just talking about
different tradeoffs since neither the call to hypot nor the respective
division is able to deliver exact results from the given arguments
unless specific conditions are met, like one of the arguments being zero
in which case we are talking about a rotation by a multiple of
90 degrees again.

Given that "clockwise" and "counterclockwise" is a lot harder to pin
down and that typical users (as you demonstrate) are likely confused
about the argument order of atan2 anyway, I don't see that this is of
the same amount of utility as specifying the transform matrix as a
rotation in degree would be.

> The hope is that it is obvious how to fix code that is using atan to
> get the angle.

If you are really using atan to get the angle, the radian version will
not actually be bad unless you use atan (or atan2) on trivial constants.

> You are right that no matter what happens, lots of code is going to
> continue to use cairo_rotate(M_PI/2), so I think fixing that function
> is a good idea.

But there is nothing to fix because the inherent problem of pi not
having a floating point representation cannot be fixed.  If you put the
assumption in that values in the vicinity of integral multiples of
M_PI/2 are supposed to be integral multiples of pi/2, you introduce
discontinuities.  Placing those discontinuities at 45+n*90 degrees
degrees is nicer than some alternatives.  But it's still messier than
providing a function working in degrees.

And if you check

git grep cairo_rotate

in the Cairo codebase, you'll find that a degree API will fix a number
of calls when it gets used.

-- 
David Kastrup