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

Submitted by David Kastrup on July 11, 2017, 6:11 p.m.

Details

Message ID 87r2xmzvw1.fsf@fencepost.gnu.org
State New
Series "Shouldn't Cairo use/offer degrees rather than radians?"
Headers show

Commit Message

David Kastrup July 11, 2017, 6:11 p.m.
Bill Spitzak <spitzak@gmail.com> writes:

> On Mon, Jul 10, 2017 at 10:10 PM, David Kastrup <dak@gnu.org> wrote:
>> Gregor Mückl <GregorMueckl@gmx.de> writes:
>>
>>> Don't pin this only on the compiler, at least on x86. The old x87
>>> compatible FPU instructions (obsolete in 64 bit mode) use internal
>>> registers with significant longer mantissa than a double. The length
>>> of the mantissa was not even fixed between implementations. However,
>>> storing values from the FPU to memory results in truncation in this
>>> case. This was "fixed" in SSE and newer intruction sets that are now
>>> fixed at 32 bits for floats and 64 bits for doubles for internal
>>> registers.
>>>
>>> I believe that this thread is starting to demonstrate that hunting for
>>> numeric precision beyond a certain point using floating point
>>> arithmetic is a fool's errand.
>>
>> The proposed code for setting angles in degrees is numerically precise
>> at right angles for any floating point representation.
>>
>> Stuff like M_PI/2 or even M_PI will always result in hunts for numeric
>> precision since pi does not have an exact floating point representation.
>>
>> This thread was about a proposed user interface that makes hunting for
>> numeric precision a non-issue for the most prominent cases while
>> representing angles in a common human-accessible format that other
>> graphic systems and formats use as well.
>>
>> So the conclusion "let's not bother offering something better because we
>> cannot conclusively figure out how bad the current situation is on
>> different systems respectively" is not really what I was pitching for.
>
> The problem is that it is not better, and providing the api will
> mislead people into thinking it is better.

This API can turn exactly by multiples of 90 degrees.  The existing API
cannot.

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

> However that is not how Cairo is implemented and changing it would
> require considerable effort for no real win. A rotate of 45 degrees is
> turned into a matrix containing sqrt(2)/2 which is not stored
> accurately.

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.

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

> Any implementation that blindly does sin(angle*(M_PI/180)) will
> produce the exact same error for right angles whether they are given
> in degrees or radians.

The patch I sent did not do this.

> So simply making the api be degrees with the easiest implementation
> will not fix anything.

DID YOU LOOK AT THE PATCH?????

To make this easier for you, I append it _again_ as an attachment to
this mail.

> I would like to make a more limited request of Cairo:
>
> 1. Change cairo_matrix_init_rotate to detect multiples of M_PI_2 and
> produce integers.

This does not work reliably because pi and consequently pi/2 do not have
exact floating point representations.  As a result, M_PI+0.5*M_PI may or
may not be equal to 1.5*M_PI.  This has to interpret values as something
that they aren't and consequently introduces errors.  You'd need to work
with tolerances, and every window of tolerance is equivalent to a window
of erroneous interpretation.  You'd need different windows for every
different floating point format.

In contrast, 0, 45, 90, ... have exact floating point representations in
any commonly used floating point format.

> From what I can see subtracting M_PI_2*n where n is
> rint(angle/M_PI_2), and then doing sin/cos works pretty good. This
> will prevent surprises in the resulting matrix. I tested some stuff in
> gnuplot and any discontinuity in the plots is below 2e-16.
>
> 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.

David Kastrup

Patch hide | download patch | download mbox

From a057b17b5018000435cc011d7a2d20072457badb Mon Sep 17 00:00:00 2001
From: David Kastrup <dak@gnu.org>
Date: Fri, 30 Jun 2017 07:01:26 +0200
Subject: [PATCH] Implement cairo_matrix_init_rotate_deg and
 cairo_matrix_rotate_deg

Those offer variants of cairo_matrix_init_rotate and cairo_matrix_rotate
that are numerically perfect at right angles without introducing
discontinuities.
---
 doc/public/cairo-sections.txt |  2 +
 src/cairo-matrix.c            | 92 +++++++++++++++++++++++++++++++++++++++++++
 src/cairo.h                   |  7 ++++
 src/cairoint.h                |  1 +
 4 files changed, 102 insertions(+)

diff --git a/doc/public/cairo-sections.txt b/doc/public/cairo-sections.txt
index 7b04ae7b3..ac5e566a0 100644
--- a/doc/public/cairo-sections.txt
+++ b/doc/public/cairo-sections.txt
@@ -430,9 +430,11 @@  cairo_matrix_init_identity
 cairo_matrix_init_translate
 cairo_matrix_init_scale
 cairo_matrix_init_rotate
+cairo_matrix_init_rotate_deg
 cairo_matrix_translate
 cairo_matrix_scale
 cairo_matrix_rotate
+cairo_matrix_rotate_deg
 cairo_matrix_invert
 cairo_matrix_multiply
 cairo_matrix_transform_distance
diff --git a/src/cairo-matrix.c b/src/cairo-matrix.c
index ae498f515..8728bcdee 100644
--- a/src/cairo-matrix.c
+++ b/src/cairo-matrix.c
@@ -312,6 +312,98 @@  cairo_matrix_rotate (cairo_matrix_t *matrix, double radians)
 }
 
 /**
+ * cairo_matrix_init_rotate_deg:
+ * @matrix: a #cairo_matrix_t
+ * @degrees: angle of rotation, in degrees. The direction of rotation
+ * is defined such that positive angles rotate in the direction from
+ * the positive X axis toward the positive Y axis. With the default
+ * axis orientation of cairo, positive angles rotate in a clockwise
+ * direction.
+ * For angles that are multiples of 90 degrees, the transformation
+ * matrix is numerically exact.
+ *
+ * Initializes @matrix to a
+ * transformation that rotates by @degrees.
+ *
+ * Since: 1.16
+ **/
+void
+cairo_matrix_init_rotate_deg (cairo_matrix_t *matrix,
+			      double degrees)
+{
+    double  s;
+    double  c;
+
+    if (unlikely (degrees <= -540.0 || degrees >= 540.0))
+      degrees = fmod (degrees, 360.0);
+    /* Now |degrees| < 540.0, and the absolute size is not larger than
+       before, so we haven't lost precision. */
+    if (unlikely (degrees <= -180.0))
+      degrees += 360.0;
+    else if (unlikely (degrees > 180.0))
+      degrees -= 360.0;
+    /* Now -180.0 < degrees <= 180.0 and we still haven't lost
+       precision.  We don't work with angles greater than 90 degrees
+       absolute in order to minimize how rounding errors of M_PI/180
+       affect the result.  The "handover" between one sine expression
+       to the next happens at angles of +-90 degrees where
+       sin(pi/2+eps) is about (1-eps^2/2).  Since the difference to
+       pi/2 should be quite small, the sine will be numerically 1 here.
+
+       Sign of the sine is chosen to avoid -0.0 in results.  This
+       version delivers exactly equal magnitude on x/y for odd
+       multiples of 45 degrees. */
+
+    if (degrees >= 0) {
+      c = sin ((90 - degrees) * (M_PI/180.0));
+      if (degrees > 90)
+	s = sin ((180 - degrees) * (M_PI/180.0));
+      else
+	s = sin (degrees * (M_PI/180.0));
+    } else {
+      c = sin ((90 + degrees) * (M_PI/180.0));
+      if (degrees < -90)
+	s = sin ((-180 - degrees) * (M_PI/180.0));
+      else
+	s = sin (degrees * (M_PI/180.0));
+    }
+
+    cairo_matrix_init (matrix,
+		       c, s,
+		       -s, c,
+		       0, 0);
+}
+slim_hidden_def(cairo_matrix_init_rotate_deg);
+
+/**
+ * cairo_matrix_rotate_deg:
+ * @matrix: a #cairo_matrix_t
+ * @degrees: angle of rotation, in degrees. The direction of rotation
+ * is defined such that positive angles rotate in the direction from
+ * the positive X axis toward the positive Y axis. With the default
+ * axis orientation of cairo, positive angles rotate in a clockwise
+ * direction.
+ * For angles that are multiples of 90 degrees, the change to the
+ * transformation matrix is numerically exact.
+ *
+ * Applies rotation by @degrees to the transformation in
+ * @matrix. The effect of the new transformation is to first rotate the
+ * coordinates by @degrees, then apply the original transformation
+ * to the coordinates.
+ *
+ * Since: 1.16
+ **/
+void
+cairo_matrix_rotate_deg (cairo_matrix_t *matrix, double degrees)
+{
+    cairo_matrix_t tmp;
+
+    cairo_matrix_init_rotate_deg (&tmp, degrees);
+
+    cairo_matrix_multiply (matrix, &tmp, matrix);
+}
+
+/**
  * cairo_matrix_multiply:
  * @result: a #cairo_matrix_t in which to store the result
  * @a: a #cairo_matrix_t
diff --git a/src/cairo.h b/src/cairo.h
index 32fc88b17..310f4c6d5 100644
--- a/src/cairo.h
+++ b/src/cairo.h
@@ -3026,6 +3026,10 @@  cairo_matrix_init_rotate (cairo_matrix_t *matrix,
 			  double radians);
 
 cairo_public void
+cairo_matrix_init_rotate_deg (cairo_matrix_t *matrix,
+			      double degrees);
+
+cairo_public void
 cairo_matrix_translate (cairo_matrix_t *matrix, double tx, double ty);
 
 cairo_public void
@@ -3034,6 +3038,9 @@  cairo_matrix_scale (cairo_matrix_t *matrix, double sx, double sy);
 cairo_public void
 cairo_matrix_rotate (cairo_matrix_t *matrix, double radians);
 
+cairo_public void
+cairo_matrix_rotate_deg (cairo_matrix_t *matrix, double degrees);
+
 cairo_public cairo_status_t
 cairo_matrix_invert (cairo_matrix_t *matrix);
 
diff --git a/src/cairoint.h b/src/cairoint.h
index 4fedf861d..f04e81131 100644
--- a/src/cairoint.h
+++ b/src/cairoint.h
@@ -1941,6 +1941,7 @@  slim_hidden_proto (cairo_mask);
 slim_hidden_proto (cairo_matrix_init);
 slim_hidden_proto (cairo_matrix_init_identity);
 slim_hidden_proto (cairo_matrix_init_rotate);
+slim_hidden_proto (cairo_matrix_init_rotate_deg);
 slim_hidden_proto (cairo_matrix_init_scale);
 slim_hidden_proto (cairo_matrix_init_translate);
 slim_hidden_proto (cairo_matrix_invert);
-- 
2.11.0


Comments

Bill Spitzak July 11, 2017, 10:42 p.m.
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) you get 1.0 and 2.220446049250313e-16 which is the same
magnitude of error.

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.

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

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


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

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. Here
it is (in python sorry):

from math import *
M_PI_2 = pi/2

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

>> 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) 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. The hope is that it is obvious how to fix
code that is using atan to get the angle. 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.
Ray Gardener July 12, 2017, 1:10 a.m.
Why not leave the Cairo API alone and implement a utility library analogous
to GLU for OpenGL and that way handy degree-based transform functions
can be put there.

Ray




On 7/11/2017 3:42 PM, Bill Spitzak wrote:
> 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) you get 1.0 and 2.220446049250313e-16 which is the same
> magnitude of error.
>
> 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.
>
>>> 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. 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.
>
>> 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.
>
>
>>> 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.
>
> 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. Here
> it is (in python sorry):
>
> from math import *
> M_PI_2 = pi/2
>
> # 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.
>
>>> 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) 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. The hope is that it is obvious how to fix
> code that is using atan to get the angle. 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.
David Kastrup July 12, 2017, 9:15 a.m.
Bill Spitzak <spitzak@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.
Bill Spitzak July 12, 2017, 9:30 p.m.
On Wed, Jul 12, 2017 at 2:15 AM, David Kastrup <dak@gnu.org> wrote:
> Bill Spitzak <spitzak@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.