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

Message ID | 87r2xmzvw1.fsf@fencepost.gnu.org |
---|---|

State | New |

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

Headers | show |

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

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.

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.

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.

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.