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

[Bill Spitzak]

On Wed, Jul 12, 2017 at 2:15 AM, David Kastrup <dak at gnu.org> wrote:
> 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)

This is all correct, however if you calculate the determinant for the matrices:

   E = sqrt(4^n/2)*eps (this is the error factor for squaring n-1 in
both examples)
   my determinant =   1 ± E^2
   your determinant = 1 ± 2E ± E^2

Your error is greater by a factor of 2E which is much larger than the
E^2 error shared by both of them.

>> 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.

Actually it appears (at least for linux glibc) that they have made a
concerted effort to make sure sin(x)^2+cos(x)^c == 1.0, apparently
sacrificing accuracy in order to achieve this (ie one sin/cos is not
actually the closest possible double representation, in order for this
identity to remain true).

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

I thought of that but it does not happen. In fact the switch happens
inside the range where it switches between returning sqrt(.5)-eps and
sqrt(.5)+eps for sin (and conversely for cos). Incrementing by epsilon
through the discontinuity and the values cleanly swap and don't go
backwards or skip any values that sin/cos can return.

>
>> 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

I get "0 0" for this.
gcc (GCC) 4.8.2 20140120 (Red Hat 4.8.2-15)

Tested with many combinations of --fast-math and -O9 and trying to
pass various numbers in argments between compilations to force
truncation to 64 bits or use of 80 bits. Although I suspect you can
get a non-zero 80 bit number, the error is less than the minimum 64
bit non-zero number.

It seems hard to believe that N*M_PI_2-N*M_PI_2 will produce anything
other than zero. If it calculated the wrong N then it would be off by
1.5 or so, not a small value like you show.

I also tried your number in my sincos function:

sincos(11*M_PI_2)
(-1.0, 0.0)

In fact I tried every integer multiple of M_PI_2 from -65556 to 65556
and did not get anything other than 1 and 0.

>> 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,

You may be right, but I certainly see a lot of programmers who are
unaware of atan2 and call atan. What I was thinking is having a
zero-brain-power translation to the correct cairo call if the
programmer is doing this.

>> 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.

You are right, the intention is for the arguments to be exactly the
same as for atan2 with y first.

> "Greater accuracy" again is a euphemism: we are just talking about

I believe the accuracy is greater than doing atan2(y,x) and then
having cairo do sin/cos. But it may be worth checking if glibc
actually gets different results for them and which one is more
correct.

> 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 intention is to give them an easy to use function if they are
*already* using atan or atan2. Though you are correct that they may be
doing cairo_rotate(-atan2(y,x)) which requires them to substitute
cairo_roate_atan2(-y,x) which is not immediately obvious.

I think you may be right that novice users are not likely to use this
to get 90 degree angles.


My conclusions:

1. Using sin/cos produces a better rotation matrix than
sin,sin(M_PI_2-x), in that the determinant is closer to 1.0. It
appears that sin/cos in the glibc libraray (and probably on the
hardware) is purposely tweaked to make this work.

2. It is nice if cairo_rotate(M_PI_2) did not put small non-zero
values in the rotation matrix as this is not expected by users and
confuses them, and may defeat some badly-written fast path detectors.
We both posted code that does this, though mine uses the better
sin/cos pair. The use of degrees is irrelevant to this, the exact same
answers are arrived at with either unit.

3. Degrees do have an advantage that some nice angles are integers
which will survive a double->float->double conversion unchanged. This
may mean that a degree api could be used to avoid some suprises. It
does not have to do anything other than convert the angle to radians
and call the radian version, the code to detect right angles can be in
the radian version. Making this division inline in the header file may
help compilers optimize and correct code that accidentally converts
radians to degrees and then calls the degree version.

4. cairo_rotate_atan2(y,x) would be a nice optimization for users
already using atan2 to generate the angle, but is not as helpful for
novices trying to get right angles as I thought.