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

Submitted by David Kastrup on June 30, 2017, 5:02 a.m.

Details

Message ID 87d19mhxub.fsf@fencepost.gnu.org
State New
Headers show
Series "Shouldn't Cairo use/offer degrees rather than radians?" ( rev: 1 ) in Cairo

Not browsing as part of any series.

Commit Message

David Kastrup June 30, 2017, 5:02 a.m.
Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:

> On Thu, 29 Jun 2017 08:08:04 +0200, David Kastrup wrote:
>
>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>> 
>>> You realize that, by introducing this distortion, you are reducing
>>> the accuracy of the computation?  
>> 
>> It replaces a sin(x)/cos(x) pair by a sin(x)/sin(pi/2-x) pair.  The
>> largest cumulative effect is at pi/4.
>
> So you introduce a kink at π/4. How is this any better, again?

There is no kink.  None at all.  The handover between different
functions happens at sin(pi/2) where there is a large stretch of 1.0.
At pi/4, there is no discontinuity but the sin/cos pair is taken as
(sin(pi/4), sin(pi/4)) instead of (sin(pi/4), cos(pi/4)).  Since pi/4
has no actual numeric representation, getting an exact angle of 45
degree in the transformation matrix may come at the cost of a larger
magnitude jitter of the determinant _iff_ you have a very, very accurate
sin/cos implementation.  Mind you: that jitter is still less than 2ulp
so you need a very good sin/cos implementation to even measure it.

Square the transform matrix for 45 degrees 20 times or so with the
degree function and the radian function and compare the results and see
which you like better.

I'll append the patch to Cairo to make it easier for you.

>>> I wonder what the term is, for deliberately producing inaccurate
>>> answers just to look good...  
>> 
>> A very well-considered tradeoff and solid engineering?
>
> “Engineering” is not the term that comes to mind. Try “marketing” or
> “management”.

Please compare the results.

> You do realize that your claim about Cairo’s “inability to reliably
> draw a half circle that seamlessly connects with another half circle
> under any rasterization” is complete nonsense, right?

Sigh.  If you cannot specify an angle for pi for which sin(angle)==0,
you will have rasterizations of half circle arcs that won't line up
perfectly because you then cannot make the half circles butt up
perfectly, and "too small to be noticeable" becomes a statistic rather
than a static feature under rasterization.

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

On Fri, 30 Jun 2017 07:02:36 +0200, David Kastrup wrote:

> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
> 
>> On Thu, 29 Jun 2017 08:08:04 +0200, David Kastrup wrote:
>>  
>>> It replaces a sin(x)/cos(x) pair by a sin(x)/sin(pi/2-x) pair.  The
>>> largest cumulative effect is at pi/4.  
>>
>> So you introduce a kink at π/4. How is this any better, again?  
> 
> There is no kink.  None at all.

No kink? Have a look--here’s your kink, at π/4.
Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:

> On Fri, 30 Jun 2017 07:02:36 +0200, David Kastrup wrote:
>
>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>> 
>>> On Thu, 29 Jun 2017 08:08:04 +0200, David Kastrup wrote:
>>>  
>>>> It replaces a sin(x)/cos(x) pair by a sin(x)/sin(pi/2-x) pair.  The
>>>> largest cumulative effect is at pi/4.  
>>>
>>> So you introduce a kink at π/4. How is this any better, again?  
>> 
>> There is no kink.  None at all.
>
> No kink? Have a look--here’s your kink, at π/4.

I have no idea what you are plotting here.  Your "fake" does not look
like a _function_ of either axis, and it also doesn't have the symmetry
of the sin/cos pairing of my code.  So if it is neither a function in
isolation nor the actual value pair, what is it?

Care to show the program you are generating your plots with?
David Kastrup <dak@gnu.org> writes:

> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>
>> On Fri, 30 Jun 2017 07:02:36 +0200, David Kastrup wrote:
>>
>>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>>> 
>>>> On Thu, 29 Jun 2017 08:08:04 +0200, David Kastrup wrote:
>>>>  
>>>>> It replaces a sin(x)/cos(x) pair by a sin(x)/sin(pi/2-x) pair.  The
>>>>> largest cumulative effect is at pi/4.  
>>>>
>>>> So you introduce a kink at π/4. How is this any better, again?  
>>> 
>>> There is no kink.  None at all.
>>
>> No kink? Have a look--here’s your kink, at π/4.
>
> I have no idea what you are plotting here.  Your "fake" does not look
> like a _function_ of either axis, and it also doesn't have the symmetry
> of the sin/cos pairing of my code.  So if it is neither a function in
> isolation nor the actual value pair, what is it?
>
> Care to show the program you are generating your plots with?

And show the units it is working with?

Cutting this a bit short: can we actually skip the part where you are
trying to show, using bravado and ad-hoc arguments, to show that I am
incompetent and don't know what I am talking about?

If you are worried about the added jitter due to the _internal_
implementation using the built-in sin/cos with radians as arguments,
I can of course add a post-correction stage (about doubling the
computational cost, of course) that will use the addition theorem in
order to get the full resolution of angles provided by the different
rasterization of degrees and radians.  But that's sort of pointless as
it will add precision amounting to something like 2^-53 .
David's code seems to be an alternative to atan2 that returns degrees.
It has a bunch of kludges to hide the fact that the atan2d he is
calling internally apparently just calls atan2 and multiplies by
180/pi. He then has another function that is in effect calculating
sin() and cos() and has even more kludges to detect multiples of 45
degrees.

The work needed in the second function indicates that passing degrees
does not really help.

Although it seems like it helps as far as I can tell if you don't want
discontinuities then adjusting the sin/cos results for N*pi/2 is
exactly the same difficulty whether degrees or radians are used.
Basically it amounts to detecting if a/90 or a/M_PI_2 is near an
integer, and is thus exactly the same problem.


On Thu, Jun 29, 2017 at 11:38 PM, David Kastrup <dak@gnu.org> wrote:
> David Kastrup <dak@gnu.org> writes:
>
>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>>
>>> On Fri, 30 Jun 2017 07:02:36 +0200, David Kastrup wrote:
>>>
>>>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>>>>
>>>>> On Thu, 29 Jun 2017 08:08:04 +0200, David Kastrup wrote:
>>>>>
>>>>>> It replaces a sin(x)/cos(x) pair by a sin(x)/sin(pi/2-x) pair.  The
>>>>>> largest cumulative effect is at pi/4.
>>>>>
>>>>> So you introduce a kink at π/4. How is this any better, again?
>>>>
>>>> There is no kink.  None at all.
>>>
>>> No kink? Have a look--here’s your kink, at π/4.
>>
>> I have no idea what you are plotting here.  Your "fake" does not look
>> like a _function_ of either axis, and it also doesn't have the symmetry
>> of the sin/cos pairing of my code.  So if it is neither a function in
>> isolation nor the actual value pair, what is it?
>>
>> Care to show the program you are generating your plots with?
>
> And show the units it is working with?
>
> Cutting this a bit short: can we actually skip the part where you are
> trying to show, using bravado and ad-hoc arguments, to show that I am
> incompetent and don't know what I am talking about?
>
> If you are worried about the added jitter due to the _internal_
> implementation using the built-in sin/cos with radians as arguments,
> I can of course add a post-correction stage (about doubling the
> computational cost, of course) that will use the addition theorem in
> order to get the full resolution of angles provided by the different
> rasterization of degrees and radians.  But that's sort of pointless as
> it will add precision amounting to something like 2^-53 .
>
> --
> David Kastrup
>
> --
> cairo mailing list
> cairo@cairographics.org
> https://lists.cairographics.org/mailman/listinfo/cairo
Attached python shows a plausible solution for making any multiple of
M_PI_2 into the expected integers. It even avoids returning negative
zero as that surprises people some times. As far as I can tell this
works for any N*M_PI_2 or (N*.5)*pi. The discontinuities near any odd
multiple of pi/4 are pretty infinitesimal as far as I can tell, and
this version is at least as accurate with magnitude as normal sin/cos.
It also has the rather nice property of guessing right angles when
extremely huge angles are put in.

However if fails if the caller accidentally does any math using float
instead of double, and this demonstrates that failure. It is possible
to fix it with the problem of adding more noise to the output by just
doing everything in float instead of double, likely this would not
produce any unwanted artifacts for Cairo users. However at very large
N, N*M_PI_2 and N*M_PI_2F get far enough apart to actually produce
different answers even when rounded to the nearest float.


On Fri, Jun 30, 2017 at 10:51 AM, Bill Spitzak <spitzak@gmail.com> wrote:
> David's code seems to be an alternative to atan2 that returns degrees.
> It has a bunch of kludges to hide the fact that the atan2d he is
> calling internally apparently just calls atan2 and multiplies by
> 180/pi. He then has another function that is in effect calculating
> sin() and cos() and has even more kludges to detect multiples of 45
> degrees.
>
> The work needed in the second function indicates that passing degrees
> does not really help.
>
> Although it seems like it helps as far as I can tell if you don't want
> discontinuities then adjusting the sin/cos results for N*pi/2 is
> exactly the same difficulty whether degrees or radians are used.
> Basically it amounts to detecting if a/90 or a/M_PI_2 is near an
> integer, and is thus exactly the same problem.
>
>
> On Thu, Jun 29, 2017 at 11:38 PM, David Kastrup <dak@gnu.org> wrote:
>> David Kastrup <dak@gnu.org> writes:
>>
>>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>>>
>>>> On Fri, 30 Jun 2017 07:02:36 +0200, David Kastrup wrote:
>>>>
>>>>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>>>>>
>>>>>> On Thu, 29 Jun 2017 08:08:04 +0200, David Kastrup wrote:
>>>>>>
>>>>>>> It replaces a sin(x)/cos(x) pair by a sin(x)/sin(pi/2-x) pair.  The
>>>>>>> largest cumulative effect is at pi/4.
>>>>>>
>>>>>> So you introduce a kink at π/4. How is this any better, again?
>>>>>
>>>>> There is no kink.  None at all.
>>>>
>>>> No kink? Have a look--here’s your kink, at π/4.
>>>
>>> I have no idea what you are plotting here.  Your "fake" does not look
>>> like a _function_ of either axis, and it also doesn't have the symmetry
>>> of the sin/cos pairing of my code.  So if it is neither a function in
>>> isolation nor the actual value pair, what is it?
>>>
>>> Care to show the program you are generating your plots with?
>>
>> And show the units it is working with?
>>
>> Cutting this a bit short: can we actually skip the part where you are
>> trying to show, using bravado and ad-hoc arguments, to show that I am
>> incompetent and don't know what I am talking about?
>>
>> If you are worried about the added jitter due to the _internal_
>> implementation using the built-in sin/cos with radians as arguments,
>> I can of course add a post-correction stage (about doubling the
>> computational cost, of course) that will use the addition theorem in
>> order to get the full resolution of angles provided by the different
>> rasterization of degrees and radians.  But that's sort of pointless as
>> it will add precision amounting to something like 2^-53 .
>>
>> --
>> David Kastrup
>>
>> --
>> cairo mailing list
>> cairo@cairographics.org
>> https://lists.cairographics.org/mailman/listinfo/cairo
Sorry I forgot to add:

David may have a point, though it is not exactly what he is stating.

Degrees have the property that "interesting" angles are more immune to
data transformations than radians. float(90) == double(90), while
float(M_PI_2) != double(M_PI_2). This can mean that you get the
expected sin and cos for "interesting" angles more often.

Where he is wrong is in thinking degrees make it easier to arrive at
these integers from an angle. I don't think it makes any difference as
you need to do an fmod or some equivalent in both cases.


On Fri, Jun 30, 2017 at 12:23 PM, Bill Spitzak <spitzak@gmail.com> wrote:
> Attached python shows a plausible solution for making any multiple of
> M_PI_2 into the expected integers. It even avoids returning negative
> zero as that surprises people some times. As far as I can tell this
> works for any N*M_PI_2 or (N*.5)*pi. The discontinuities near any odd
> multiple of pi/4 are pretty infinitesimal as far as I can tell, and
> this version is at least as accurate with magnitude as normal sin/cos.
> It also has the rather nice property of guessing right angles when
> extremely huge angles are put in.
>
> However if fails if the caller accidentally does any math using float
> instead of double, and this demonstrates that failure. It is possible
> to fix it with the problem of adding more noise to the output by just
> doing everything in float instead of double, likely this would not
> produce any unwanted artifacts for Cairo users. However at very large
> N, N*M_PI_2 and N*M_PI_2F get far enough apart to actually produce
> different answers even when rounded to the nearest float.
>
>
> On Fri, Jun 30, 2017 at 10:51 AM, Bill Spitzak <spitzak@gmail.com> wrote:
>> David's code seems to be an alternative to atan2 that returns degrees.
>> It has a bunch of kludges to hide the fact that the atan2d he is
>> calling internally apparently just calls atan2 and multiplies by
>> 180/pi. He then has another function that is in effect calculating
>> sin() and cos() and has even more kludges to detect multiples of 45
>> degrees.
>>
>> The work needed in the second function indicates that passing degrees
>> does not really help.
>>
>> Although it seems like it helps as far as I can tell if you don't want
>> discontinuities then adjusting the sin/cos results for N*pi/2 is
>> exactly the same difficulty whether degrees or radians are used.
>> Basically it amounts to detecting if a/90 or a/M_PI_2 is near an
>> integer, and is thus exactly the same problem.
>>
>>
>> On Thu, Jun 29, 2017 at 11:38 PM, David Kastrup <dak@gnu.org> wrote:
>>> David Kastrup <dak@gnu.org> writes:
>>>
>>>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>>>>
>>>>> On Fri, 30 Jun 2017 07:02:36 +0200, David Kastrup wrote:
>>>>>
>>>>>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>>>>>>
>>>>>>> On Thu, 29 Jun 2017 08:08:04 +0200, David Kastrup wrote:
>>>>>>>
>>>>>>>> It replaces a sin(x)/cos(x) pair by a sin(x)/sin(pi/2-x) pair.  The
>>>>>>>> largest cumulative effect is at pi/4.
>>>>>>>
>>>>>>> So you introduce a kink at π/4. How is this any better, again?
>>>>>>
>>>>>> There is no kink.  None at all.
>>>>>
>>>>> No kink? Have a look--here’s your kink, at π/4.
>>>>
>>>> I have no idea what you are plotting here.  Your "fake" does not look
>>>> like a _function_ of either axis, and it also doesn't have the symmetry
>>>> of the sin/cos pairing of my code.  So if it is neither a function in
>>>> isolation nor the actual value pair, what is it?
>>>>
>>>> Care to show the program you are generating your plots with?
>>>
>>> And show the units it is working with?
>>>
>>> Cutting this a bit short: can we actually skip the part where you are
>>> trying to show, using bravado and ad-hoc arguments, to show that I am
>>> incompetent and don't know what I am talking about?
>>>
>>> If you are worried about the added jitter due to the _internal_
>>> implementation using the built-in sin/cos with radians as arguments,
>>> I can of course add a post-correction stage (about doubling the
>>> computational cost, of course) that will use the addition theorem in
>>> order to get the full resolution of angles provided by the different
>>> rasterization of degrees and radians.  But that's sort of pointless as
>>> it will add precision amounting to something like 2^-53 .
>>>
>>> --
>>> David Kastrup
>>>
>>> --
>>> cairo mailing list
>>> cairo@cairographics.org
>>> https://lists.cairographics.org/mailman/listinfo/cairo
A function like this might be useful:

  cairo_rotate_2(cairo_t*, double x, double y)

This would rotate it so the x axis is angled by atan2(y,x). I am not
sure if it should also scale by hypot(x,y) or normalize the vector.
This can be implemented without any trig functions.


On Fri, Jun 30, 2017 at 12:31 PM, Bill Spitzak <spitzak@gmail.com> wrote:
> Sorry I forgot to add:
>
> David may have a point, though it is not exactly what he is stating.
>
> Degrees have the property that "interesting" angles are more immune to
> data transformations than radians. float(90) == double(90), while
> float(M_PI_2) != double(M_PI_2). This can mean that you get the
> expected sin and cos for "interesting" angles more often.
>
> Where he is wrong is in thinking degrees make it easier to arrive at
> these integers from an angle. I don't think it makes any difference as
> you need to do an fmod or some equivalent in both cases.
>
>
> On Fri, Jun 30, 2017 at 12:23 PM, Bill Spitzak <spitzak@gmail.com> wrote:
>> Attached python shows a plausible solution for making any multiple of
>> M_PI_2 into the expected integers. It even avoids returning negative
>> zero as that surprises people some times. As far as I can tell this
>> works for any N*M_PI_2 or (N*.5)*pi. The discontinuities near any odd
>> multiple of pi/4 are pretty infinitesimal as far as I can tell, and
>> this version is at least as accurate with magnitude as normal sin/cos.
>> It also has the rather nice property of guessing right angles when
>> extremely huge angles are put in.
>>
>> However if fails if the caller accidentally does any math using float
>> instead of double, and this demonstrates that failure. It is possible
>> to fix it with the problem of adding more noise to the output by just
>> doing everything in float instead of double, likely this would not
>> produce any unwanted artifacts for Cairo users. However at very large
>> N, N*M_PI_2 and N*M_PI_2F get far enough apart to actually produce
>> different answers even when rounded to the nearest float.
>>
>>
>> On Fri, Jun 30, 2017 at 10:51 AM, Bill Spitzak <spitzak@gmail.com> wrote:
>>> David's code seems to be an alternative to atan2 that returns degrees.
>>> It has a bunch of kludges to hide the fact that the atan2d he is
>>> calling internally apparently just calls atan2 and multiplies by
>>> 180/pi. He then has another function that is in effect calculating
>>> sin() and cos() and has even more kludges to detect multiples of 45
>>> degrees.
>>>
>>> The work needed in the second function indicates that passing degrees
>>> does not really help.
>>>
>>> Although it seems like it helps as far as I can tell if you don't want
>>> discontinuities then adjusting the sin/cos results for N*pi/2 is
>>> exactly the same difficulty whether degrees or radians are used.
>>> Basically it amounts to detecting if a/90 or a/M_PI_2 is near an
>>> integer, and is thus exactly the same problem.
>>>
>>>
>>> On Thu, Jun 29, 2017 at 11:38 PM, David Kastrup <dak@gnu.org> wrote:
>>>> David Kastrup <dak@gnu.org> writes:
>>>>
>>>>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>>>>>
>>>>>> On Fri, 30 Jun 2017 07:02:36 +0200, David Kastrup wrote:
>>>>>>
>>>>>>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>>>>>>>
>>>>>>>> On Thu, 29 Jun 2017 08:08:04 +0200, David Kastrup wrote:
>>>>>>>>
>>>>>>>>> It replaces a sin(x)/cos(x) pair by a sin(x)/sin(pi/2-x) pair.  The
>>>>>>>>> largest cumulative effect is at pi/4.
>>>>>>>>
>>>>>>>> So you introduce a kink at π/4. How is this any better, again?
>>>>>>>
>>>>>>> There is no kink.  None at all.
>>>>>>
>>>>>> No kink? Have a look--here’s your kink, at π/4.
>>>>>
>>>>> I have no idea what you are plotting here.  Your "fake" does not look
>>>>> like a _function_ of either axis, and it also doesn't have the symmetry
>>>>> of the sin/cos pairing of my code.  So if it is neither a function in
>>>>> isolation nor the actual value pair, what is it?
>>>>>
>>>>> Care to show the program you are generating your plots with?
>>>>
>>>> And show the units it is working with?
>>>>
>>>> Cutting this a bit short: can we actually skip the part where you are
>>>> trying to show, using bravado and ad-hoc arguments, to show that I am
>>>> incompetent and don't know what I am talking about?
>>>>
>>>> If you are worried about the added jitter due to the _internal_
>>>> implementation using the built-in sin/cos with radians as arguments,
>>>> I can of course add a post-correction stage (about doubling the
>>>> computational cost, of course) that will use the addition theorem in
>>>> order to get the full resolution of angles provided by the different
>>>> rasterization of degrees and radians.  But that's sort of pointless as
>>>> it will add precision amounting to something like 2^-53 .
>>>>
>>>> --
>>>> David Kastrup
>>>>
>>>> --
>>>> cairo mailing list
>>>> cairo@cairographics.org
>>>> https://lists.cairographics.org/mailman/listinfo/cairo
Bill Spitzak <spitzak@gmail.com> writes:

> David's code seems to be an alternative to atan2 that returns degrees.

The patch I posted does not do anything akin to atan2.

> It has a bunch of kludges to hide the fact that the atan2d he is
> calling internally

I am not calling any atan2d.

> apparently just calls atan2 and multiplies by 180/pi. He then has
> another function that is in effect calculating sin() and cos() and has
> even more kludges to detect multiples of 45 degrees.

Are you kidding me?  The code I posted does not do anything with
multiples of 45 degrees.  It does not detect anything either: it figures
out the quadrant the angle is in while still in degrees and then uses
the pi/180 factor only with respect to the given quadrant in order to
minimize the accumulative error caused by pi/180 not having an exact
representation to start with.

> The work needed in the second function indicates that passing degrees
> does not really help.

A few comparisons?

> Although it seems like it helps as far as I can tell if you don't want
> discontinuities then adjusting the sin/cos results for N*pi/2 is
> exactly the same difficulty whether degrees or radians are used.

N*pi/2 does not have an exact representation.  3*pi/2 may or may not be
pi+pi/2.  pi+pi/2 may or may not be 3*pi/2.  pi+pi/2+pi/2 is unlikely to
be 2*pi.  That is _not_, I repeat _not_ "exactly the same difficulty" as
working with multiples of 90 which are perfectly representable in all
floating point formats in use (9 bits of mantissa will suffice for
binary floating point, and decimal floating point works too).

> Basically it amounts to detecting if a/90 or a/M_PI_2 is near an
> integer, and is thus exactly the same problem.

Sorry, but that's just wrong.  Do the same code I did for degrees for
radians, and it just won't work.  You won't get dependable zeros for
right angles because right angles don't have a floating point
representation to start with in radians.
Bill Spitzak <spitzak@gmail.com> writes:

> Sorry I forgot to add:
>
> David may have a point, though it is not exactly what he is stating.
>
> Degrees have the property that "interesting" angles are more immune to
> data transformations than radians. float(90) == double(90), while
> float(M_PI_2) != double(M_PI_2).

The problem rather is that any M_PI_2 is not equal to pi/2.  There is no
floating point number p2 for which cos(p2) == 0.0.

> This can mean that you get the expected sin and cos for "interesting"
> angles more often.

Well, more often than _zero_ times is not a high bar, is it?

> Where he is wrong is in thinking degrees make it easier to arrive at
> these integers from an angle. I don't think it makes any difference as
> you need to do an fmod or some equivalent in both cases.

So for what value of p is cos(p) == 0.0 ?
Bill Spitzak <spitzak@gmail.com> writes:

> A function like this might be useful:
>
>   cairo_rotate_2(cairo_t*, double x, double y)
>
> This would rotate it so the x axis is angled by atan2(y,x). I am not
> sure if it should also scale by hypot(x,y) or normalize the vector.
> This can be implemented without any trig functions.

Not much more user-friendly than setting the transform matrix to
individual values.
On Fri, Jun 30, 2017 at 1:44 PM, David Kastrup <dak@gnu.org> wrote:

> I am not calling any atan2d.

Your code copied and pasted from the previous email. If this is not
your code please explain where the code you want to use is:

Real
Offset::angle_degrees () const
{
  Real x = coordinate_a_ [X_AXIS];
  Real y = coordinate_a_ [Y_AXIS];

  // We keep in the vicinity of multiples of 45 degrees here: this is
  // where straightforward angles for straightforward angular
  // relations are most expected.  The factors of 2 employed in the
  // comparison are not really perfect for that: sqrt(2)+1 would be
  // the factor giving exact windows of 45 degrees rather than what we
  // have here.  It's just that 2 is likely to generate nicer code
  // than 2.4 and the exact handover does not really matter.
  //
  // Comparisons here are chosen appropriately to let infinities end
  // up in their "exact" branch.  As opposed to the normal atan2
  // function behavior, this makes "competing" infinities result in
  // NAN angles.
  if (y < 0.0)
    {
      if (2*x < -y)
        if (-x > -2*y)          // x < 0, y < 0, |x| > |2y|
          return -180 + atan2d (-y, -x);
        else if (-2*x >= -y)    // x < 0, y < 0, |y| < |2x| <= |4y|
          return -135 + atan2d (x - y, -y - x);
        else                    // y < 0, |y| >= |2x|
          return -90 + atan2d (x, -y);
      else if (x <= -2*y)       // x > 0, y < 0, |y| <= |2x| < |4y|
        return -45 + atan2d (x + y, x - y);
      // Drop through for y < 0, x > |2y|
    }
  else if (y > 0.0)
    {
      if (2*x < y)
        if (-x > 2*y)           // x < 0, y >= 0, |x| > |2y|
          return 180 - atan2d (y, -x);
        else if (-2*x >= y)     // x < 0, y >= 0, |y| < |2x| <= |4y|
          return 135 - atan2d (x + y, y - x);
        else                    // y >= 0, |y| >= |2x|
          return 90 - atan2d (x, y);
      else if (x <= 2*y)        // x >= 0, y >= 0, |y| < |2x| < |4y|
        return 45 - atan2d (x - y, x + y);
      // Drop through for y > 0, x > |2y|
    }
  else
    // we return 0 for (0,0).  NAN would be an option but is a
    // nuisance for getting back to rectangular coordinates.  Strictly
    // speaking, this argument would be just as valid for (+inf.0,
    // +inf.0), but then infinities are already an indication of a
    // problem in LilyPond.
    return (x < 0.0) ? 180 : 0;
  return atan2d (y, x);
}
On Fri, 30 Jun 2017 07:02:36 +0200, David Kastrup wrote:

> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
> 
>> You do realize that your claim about Cairo’s “inability to reliably
>> draw a half circle that seamlessly connects with another half circle
>> under any rasterization” is complete nonsense, right?  
> 
> Sigh.  If you cannot specify an angle for pi for which sin(angle)==0,
> you will have rasterizations of half circle arcs that won't line up
> perfectly because you then cannot make the half circles butt up
> perfectly, and "too small to be noticeable" becomes a statistic rather
> than a static feature under rasterization.

There is no such situation.
Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:

> On Fri, 30 Jun 2017 07:02:36 +0200, David Kastrup wrote:
>
>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>> 
>>> You do realize that your claim about Cairo’s “inability to reliably
>>> draw a half circle that seamlessly connects with another half circle
>>> under any rasterization” is complete nonsense, right?  
>> 
>> Sigh.  If you cannot specify an angle for pi for which sin(angle)==0,
>> you will have rasterizations of half circle arcs that won't line up
>> perfectly because you then cannot make the half circles butt up
>> perfectly, and "too small to be noticeable" becomes a statistic rather
>> than a static feature under rasterization.
>
> There is no such situation.
What is wrong with this list?
2017-07-03 17:14 GMT+02:00 David Kastrup <dak@gnu.org>:
> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>
>> On Fri, 30 Jun 2017 07:02:36 +0200, David Kastrup wrote:
>>
>>> Lawrence D'Oliveiro <ldo@geek-central.gen.nz> writes:
>>>
>>>> You do realize that your claim about Cairo’s “inability to reliably
>>>> draw a half circle that seamlessly connects with another half circle
>>>> under any rasterization” is complete nonsense, right?
>>>
>>> Sigh.  If you cannot specify an angle for pi for which sin(angle)==0,
>>> you will have rasterizations of half circle arcs that won't line up
>>> perfectly because you then cannot make the half circles butt up
>>> perfectly, and "too small to be noticeable" becomes a statistic rather
>>> than a static feature under rasterization.
>>
>> There is no such situation.
>
> What is wrong with this list?

I don't think it is the list as a whole. But I do understand your
frustration at this point.

Guillermo
Gregor Mückl <GregorMueckl@gmx.de> writes:

> On 07/03/17 17:14, David Kastrup wrote:
>> [example code for drawing a half-circle omitted]
>
> Could you please explain why you chose that y value of
> 0.501954078674316 specifically? At first sight, it seems to be
> specifically designed to hit a pathological case.

Of course, because of the particular preceding discussion with Lawrence
D'Oliveiro.  I quote:

    >> Sigh.  If you cannot specify an angle for pi for which
    >> sin(angle)==0, you will have rasterizations of half circle arcs
    >> that won't line up perfectly because you then cannot make the
    >> half circles butt up perfectly, and "too small to be noticeable"
    >> becomes a statistic rather than a static feature under
    >> rasterization.
    >
    > There is no such situation.

This is a counterexample to the "there is no such situation claim" of
Lawrence.

> How big of an epsilon can you add to/subtract from this value before
> the output changes?

Quite small, of course: this value was bisected until failure.  I knew
that failure was going to happen (just like Lawrence knew that it
couldn't), so I took a half circle and chased the edge case down.  The
window is small since its size correlates to the inability of "double"
to properly represent pi.  It will be different for different floating
point implementations/representations.

The point is that radians are unable to represent important angles for
graphics operations in normal floating point formats and that there are
consequences.  They are more obvious for transform matrices but they are
even there for arcs.

> I haven't played around with the program.

It's an edge case.  Not more, but also not less.
On 07/03/17 17:14, David Kastrup wrote:
 > [example code for drawing a half-circle omitted]

Could you please explain why you chose that y value of 0.501954078674316 
specifically? At first sight, it seems to be specifically designed to 
hit a pathological case. How big of an epsilon can you add to/subtract 
from this value before the output changes? I haven't played around with 
the program.

Gregor
Hello colleague(s),

i actually wanted to look closer on the example, unfortunately i
cannot reproduce. This is a linux box, intel i3, ubuntu 14.04.

Can you explain further, how you decided for the y- coordinate?
Andreas Lobinger <lobingera@gmail.com> writes:

> Hello colleague(s),
>
> i actually wanted to look closer on the example, unfortunately i
> cannot reproduce. This is a linux box, intel i3, ubuntu 14.04.

Intel Penryn (Core Duo) and Ubuntu 17.10 (to be) here.

> Can you explain further, how you decided for the y- coordinate?

Bisection: make the value smaller when the half circle did not touch the
border (like in your picture), make it larger when it did touch.

Eventually I reached a value where only one leg touched.  14.04 is a
long time ago: conceivable that there were changes in rasterization
since then.  Still I consider it rather likely that you'll find such a
half-touching value.

I did not try to automate this.
Hello colleagues(s),

i checked, i run a 1.13.1, but 1.14.6 gets me the same result.

btw: if i run this arc through cairo_copy_path and readout the path
data (see example in the API), i get:

moveto[20.0,0.503906]
curveto[20.0,7.16797,10.0,7.16797,10.0,0.503906]

rounding of digits is a side effect, the y coordinates on the moveto
and curveto are the same.
Andreas Lobinger <lobingera@gmail.com> writes:

> Hello colleagues(s),
>
> i checked, i run a 1.13.1, but 1.14.6 gets me the same result.
>
> btw: if i run this arc through cairo_copy_path and readout the path
> data (see example in the API), i get:
>
> moveto[20.0,0.503906]
> curveto[20.0,7.16797,10.0,7.16797,10.0,0.503906]
>
> rounding of digits is a side effect, the y coordinates on the moveto
> and curveto are the same.

With this print format, the numbers would look the same.  Do they
compare equal when comparing them with == ?  And/or can you print with
something like 20 significant digits?
David Kastrup <dak@gnu.org> writes:

> Andreas Lobinger <lobingera@gmail.com> writes:
>
>> Hello colleagues(s),
>>
>> i checked, i run a 1.13.1, but 1.14.6 gets me the same result.
>>
>> btw: if i run this arc through cairo_copy_path and readout the path
>> data (see example in the API), i get:
>>
>> moveto[20.0,0.503906]
>> curveto[20.0,7.16797,10.0,7.16797,10.0,0.503906]
>>
>> rounding of digits is a side effect, the y coordinates on the moveto
>> and curveto are the same.
>
> With this print format, the numbers would look the same.  Do they
> compare equal when comparing them with == ?  And/or can you print with
> something like 20 significant digits?

The reason I ask: if the y coordinates of the curveto are _exactly_
identical, this artifact would not actually be a consequence of
trigonometric functions unless cairo_copy_path is responsible for
minuscule changes as well.
Hello colleague(s),

On Wed, Jul 5, 2017 at 10:37 PM, David Kastrup <dak@gnu.org> wrote:

> With this print format, the numbers would look the same.  Do they
> compare equal when comparing them with == ?  And/or can you print with
> something like 20 significant digits?
>
> The reason I ask: if the y coordinates of the curveto are _exactly_
> identical, this artifact would not actually be a consequence of
> trigonometric functions unless cairo_copy_path is responsible for
> minuscule changes as well.

step1: i'm attaching the c-code i used to print the data
step2: output (with %.20f)

lobi@orange4:~/juliarepo$ gcc attachment1.c -I/usr/include/cairo/ -lcairo
lobi@orange4:~/juliarepo$ ./a.out
moveto  20.00000000000000000000,0.50390625000000000000
curveto 20.00000000000000000000,7.16796875000000000000
10.00000000000000000000,7.16796875000000000000
10.00000000000000000000,0.50390625000000000000

(and i'm pretty much convinced that path_copy just copies data without
interfering)

Without further debugging or looking into the code i'd assume we see
here the process of cairo modeling the arc piecewise by curveto
operators. And in this case a single curve - regarding the scaling and
resolution of the output device seems to be enough.

step3: please run and post. I'm really interested how your cairo is
creating a different curve.
I had to change the format to CAIRO_FORMAT_ARGB32 to get this to work
(otherwise I got CAIRO_STATUS_INVALID_FORMAT) and my output image is
"correct" in that both ends are the same y pixel. The path prints
exactly the same numbers as you show.

I suspect fiddling with your y value would get the error back however.

On Thu, Jul 6, 2017 at 10:37 AM, Andreas Lobinger <lobingera@gmail.com> wrote:
> Hello colleague(s),
>
> On Wed, Jul 5, 2017 at 10:37 PM, David Kastrup <dak@gnu.org> wrote:
>
>> With this print format, the numbers would look the same.  Do they
>> compare equal when comparing them with == ?  And/or can you print with
>> something like 20 significant digits?
>>
>> The reason I ask: if the y coordinates of the curveto are _exactly_
>> identical, this artifact would not actually be a consequence of
>> trigonometric functions unless cairo_copy_path is responsible for
>> minuscule changes as well.
>
> step1: i'm attaching the c-code i used to print the data
> step2: output (with %.20f)
>
> lobi@orange4:~/juliarepo$ gcc attachment1.c -I/usr/include/cairo/ -lcairo
> lobi@orange4:~/juliarepo$ ./a.out
> moveto  20.00000000000000000000,0.50390625000000000000
> curveto 20.00000000000000000000,7.16796875000000000000
> 10.00000000000000000000,7.16796875000000000000
> 10.00000000000000000000,0.50390625000000000000
>
> (and i'm pretty much convinced that path_copy just copies data without
> interfering)
>
> Without further debugging or looking into the code i'd assume we see
> here the process of cairo modeling the arc piecewise by curveto
> operators. And in this case a single curve - regarding the scaling and
> resolution of the output device seems to be enough.
>
> step3: please run and post. I'm really interested how your cairo is
> creating a different curve.
>
> --
> cairo mailing list
> cairo@cairographics.org
> https://lists.cairographics.org/mailman/listinfo/cairo
If I take a look at the pixels, I can see the arc not making a half circle in this case. When I run the following program I get this output. This is with cairo 1.14.6. If I try the same program with Valgrind the half circle lines up?

Eric

...
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0 255 255 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
255 255 255 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
255 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 

...

With Valgrind
...
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0 255 255 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0 255 255 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0 255 255   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
...

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

int main()
{
  //cairo_surface_t *surface = cairo_image_surface_create (CAIRO_FORMAT_RGB16_565, 30, 20);
  cairo_surface_t *surface = cairo_image_surface_create (CAIRO_FORMAT_ARGB32, 30, 20);
  cairo_t *cr = cairo_create (surface);
  cairo_set_line_width (cr, 2);
  cairo_set_antialias (cr, CAIRO_ANTIALIAS_NONE);
  cairo_set_source_rgb (cr, 0.0, 0.0, 0.0);
  cairo_paint(cr);

  cairo_set_source_rgb (cr, 0.0, 1.0, 0.0);
  cairo_arc (cr, 15, 0.501954078674316, 5, 0, M_PI);
  cairo_stroke (cr);
  cairo_destroy (cr);

  int i=0;
  int j=0;
  cairo_surface_flush(surface);
  int width = cairo_image_surface_get_width(surface);
  int height = cairo_image_surface_get_height(surface);
  int stride = cairo_format_stride_for_width (CAIRO_FORMAT_ARGB32, 30);
  int channels = 4;
  unsigned char *s=cairo_image_surface_get_data(surface);
  unsigned char *pixels=s;
  for(i=0;i<width;i++)
    {
      for(j=0;j<height;j++)
        {
          printf("%3i ", pixels[j*stride+i*channels+1]);
        }
      printf("\n");
    }

  cairo_surface_destroy (surface);
  return 0;
}
Has anyone figured out a good explanation for this? I haven't been able to figure it out. It looks like the problem is tied to only one particular precision point. I don't know why that is. If I change the size of the surface, radius, y-point or even pi I get the same problem point. If I run with valgrind I get a different repetitive problem point. Maybe on some hardware this doesn't happen? I am testing on a Intel 32 bit Atom CPU N270 1.60GHz × 2 with Ubuntu16.04. This is my latest to try to figure out what is happening. 

Eric


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

int main()
{
  printf("%s\n", cairo_version_string());
  cairo_surface_t *surface = cairo_image_surface_create (CAIRO_FORMAT_ARGB32, 30, 20);
  cairo_t *cr = cairo_create (surface);
  cairo_set_line_width (cr, 2);
  cairo_set_antialias (cr, CAIRO_ANTIALIAS_NONE);

  double test_number=0.0;
  double increment=0.1;
  int channels=4;
  int width=cairo_image_surface_get_width(surface);
  int height=cairo_image_surface_get_height(surface);
  int stride=cairo_format_stride_for_width (CAIRO_FORMAT_ARGB32, width);
  printf("width %i, height %i, stride %i\n", width, height, stride);
  unsigned char *s=cairo_image_surface_get_data(surface);
  unsigned char *pixels=s;
  //Safety for do loop.
  int safety=0;
  int counter=0;
  int i=0;
  int j=0;
  
  //Change the arc radius for each loop. Use i to change the test number for new search also.
  for(i=0;i<5;i++)
    {
      do
        {
          test_number+=increment;      
          cairo_set_source_rgb(cr, 0.0, 0.0, 0.0);
          cairo_paint(cr);
          cairo_set_source_rgb(cr, 0.0, 1.0, 0.0);
          //Test changing arc radius to see if this changes the numbers.
          cairo_arc(cr, 15, test_number, 5+i, 0, M_PI);
          cairo_stroke(cr);  
          cairo_surface_flush(surface); 

          //Look for the color pixels in the row.
          counter=0;
          for(j=0;j<width;j++)
            {
              if(pixels[i*stride+j*channels+1]!=0) counter++;
            } 
          //Refine precision if the row was all 0's. Backup and start looking again.
          if(counter==0)
            {
              //Go back one.
              test_number-=increment;
              increment=increment*0.1;
              printf("%i %i %i %.15f, %.15f\n", i, counter, safety, test_number, increment);
            }
         
          if(safety>100)
            {
              printf("Couldn't find partial arc.\n");
              break;
            }
          safety++;    
        }while(counter<1||counter>3);

      printf("Row %i %.15f ", i, test_number);
      for(j=0;j<width;j++)
        {
          printf("%i, ", pixels[i*stride+j*channels+1]);
        }
      printf("\n");

      safety=0;
      test_number=(double)i+1;
      increment=0.1;
    }

  cairo_destroy(cr);
  cairo_surface_destroy(surface);
  cairo_debug_reset_static_data();
  return 0;
}
cecashon@aol.com writes:

> Has anyone figured out a good explanation for this? I haven't been
> able to figure it out. It looks like the problem is tied to only one
> particular precision point. I don't know why that is. If I change the
> size of the surface, radius, y-point or even pi I get the same problem
> point. If I run with valgrind I get a different repetitive problem
> point.

Does valgrind link debug libraries in some manner?

> Maybe on some hardware this doesn't happen?

That's certainly quite likely since it depends on just what sin/cos do
on a minuscule scale where different FPUs will quite likely deliver
different results.

> I am testing on a Intel 32 bit Atom CPU N270 1.60GHz × 2 with
> Ubuntu16.04. This is my latest to try to figure out what is happening.

Compilers may choose to use more floating point precision than specified
for intermediate results.  That means that different GCC versions and
different optimization options might change values around.

For that reason it's important that you don't link with any different
libraries, even libraries compiled from the same source code.
One can do cairo_matrix_init right? Just directly set unitary values
to get perfect right-angle rotations?

Ray


On 7/10/2017 1:04 PM, David Kastrup wrote:
> cecashon@aol.com writes:
>
>> Has anyone figured out a good explanation for this? I haven't been
>> able to figure it out. It looks like the problem is tied to only one
>> particular precision point. I don't know why that is. If I change the
>> size of the surface, radius, y-point or even pi I get the same problem
>> point. If I run with valgrind I get a different repetitive problem
>> point.
> Does valgrind link debug libraries in some manner?
>
>> Maybe on some hardware this doesn't happen?
> That's certainly quite likely since it depends on just what sin/cos do
> on a minuscule scale where different FPUs will quite likely deliver
> different results.
>
>> I am testing on a Intel 32 bit Atom CPU N270 1.60GHz × 2 with
>> Ubuntu16.04. This is my latest to try to figure out what is happening.
> Compilers may choose to use more floating point precision than specified
> for intermediate results.  That means that different GCC versions and
> different optimization options might change values around.
>
> For that reason it's important that you don't link with any different
> libraries, even libraries compiled from the same source code.
>
On 07/10/2017 10:04 PM, David Kastrup wrote:
> cecashon@aol.com writes:
> 
>> Has anyone figured out a good explanation for this? I haven't been
>> able to figure it out. It looks like the problem is tied to only one
>> particular precision point. I don't know why that is. If I change the
>> size of the surface, radius, y-point or even pi I get the same problem
>> point. If I run with valgrind I get a different repetitive problem
>> point.
> 
> Does valgrind link debug libraries in some manner?
> 

No, valgrind emulates x86 completely in software. It does not load 
different libraries, but may hook certain libc function calls to perform 
its job.

>> Maybe on some hardware this doesn't happen?
> 
> That's certainly quite likely since it depends on just what sin/cos do
> on a minuscule scale where different FPUs will quite likely deliver
> different results.
> 
>> I am testing on a Intel 32 bit Atom CPU N270 1.60GHz × 2 with
>> Ubuntu16.04. This is my latest to try to figure out what is happening.
> 
> Compilers may choose to use more floating point precision than specified
> for intermediate results.  That means that different GCC versions and
> different optimization options might change values around.
> 
> For that reason it's important that you don't link with any different
> libraries, even libraries compiled from the same source code.
> 

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.
Ray Gardener <rayg@daylongraphics.com> writes:

> One can do cairo_matrix_init right? Just directly set unitary values
> to get perfect right-angle rotations?

You can see in the Cairo code base itself (admittedly test and perf
directories rather than the core but they are typical basic uses) that
people prefer not to even when it would help precision.
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.
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.

You wanted two rotates by 45 degrees to be perfect. 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. The
exact same value is stored whether the rotation is sent as degrees or
radians.

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. So simply making the api be degrees with the
easiest implementation will not fix anything.

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