modify coord-rotate to get exact values for (sin PI) etc

modify coord-rotate to get exact values for (sin PI) etc

issue 4640

Done by switching to appropiate values for the angle and/or switching
sin to cos and vice versa

[thomasmorley651, 2015-10-19 23:15:04]

please review

[dak, 2015-10-20 14:25:21]

I'm not really enthused about this approach even though it originates from a
proposal of mine.  It works basically by moving around problematic angles,
leading to discontinuities at multiples of PI/4 rather than just PI being
imprecise.

Fundamentally, we cannot fix this properly at the Scheme layer: even in
GUILE-2.0, we get

scheme@(guile-user)> (sin (atan 0 -1))
$1 = 1.2246467991473532e-16

I think what we should be aiming for is grounding LilyPond's offering of
trigonometric operators in degrees.

This is what PostScript offers:

GS<2>dak@lola:/usr/local/tmp/lilypond/elisp$ gsnd
GPL Ghostscript 9.15 (2014-09-22)
Copyright (C) 2014 Artifex Software, Inc.  All rights reserved.
This software comes with NO WARRANTY: see the file PUBLIC for details.
GS>-1 0 atan ==
270.0
GS>0 -1 atan ==
180.0
GS>0 -1 atan sin ==
0.0
GS>quit

This is also how METAFONT rolls:
dak@lola:/usr/local/tmp/lilypond/elisp$ mf '\\relax'
This is METAFONT, Version 2.7182818 (TeX Live 2015/dev/Debian) (preloaded
base=mf)

*show angle(-1,0)

*;
>> 180
*show sin(angle(-1,0));
>> sin
! Extra tokens will be flushed.
<to be read again> 
                   (
<*> show sin(
             angle(-1,0));
? 

*show sin;
>> sin
*show sind(angle(-1,0));
>> 0
*

If those graphical languages that we also work with decide it to be best to keep
off radians, maybe we should take a hint from them.

Our degree engine might even end up quite similar to what we have here, just
working for multiples of 45 rather than PI/4, whether implemented in Scheme or
C.  But the major difference is that the computer is _guaranteed_ to have an
exact representation of multiples of 45.  In contrast, we have no guarantee that
3*(PI/4) has an exact representation: this will only be the case if the lowest
two significant bits of the representation of PI/4 are zero since 4 > 3pi/4 > 2
and 1 > pi/4 > 0.5 and thus the mantissa of 3pi/4 takes up two more bits in the
most significant bits.

So we can't hope to keep stuff exact when juggling with multiples of angles in
radians.

So let's rather discuss which degree-based primitives we'd rather have,
implement them, and replace everything in LilyPond juggling with angles with
them and make sure that our own commands only take degrees rather than radians.

But I fear that this radian-based approach is only half a bandaid that will
result in this problem reappearing in different forms.

[thomasmorley651, 2015-10-20 21:55:00]

Let me split my reply into two.

This one about the current patch and how to proceed.

On 2015/10/20 14:25:21, dak wrote:
> I'm not really enthused about this approach even though it originates from a
> proposal of mine.  It works basically by moving around problematic angles,

exactly

> leading to discontinuities at multiples of PI/4 rather than just PI being
> imprecise.

PI/2 is problematic as well (from scheme-sandbox):
guile> (cos (/ PI 2))
6.12323399573677e-17


I locally have a bunch of code, which I plan to turn into patches, committing
them in small pieces to ease review.

- simplified make-bezier-sandwich-stencil (already pushed)
- exact coord-rotate 
- a generic make-bow-stencil, which _needs_ _exact_ coord-rotate
- the undertie-markup-command
- redefined make-parentheses-stencil
- a markup-command for drawing squiggle lines (under construction)
- a stencil-command to draw more fancy squiggle lines (under construction)

All rely on make-bow-stencil.


I could 
(a)
upload them, one after another, postponing a better coord-rotate (but again, an
exact coord-rotate for multiples of PI/2 is needed)
or 
(b)
doing a better coord-rotate first (see next post)

What do you think?

[thomasmorley651, 2015-10-20 22:10:50]

This oe is about radians vs degrees

On 2015/10/20 14:25:21, dak wrote:
> Fundamentally, we cannot fix this properly at the Scheme layer: even in
> GUILE-2.0, we get
> 
> scheme@(guile-user)> (sin (atan 0 -1))
> $1 = 1.2246467991473532e-16
> 
> I think what we should be aiming for is grounding LilyPond's offering of
> trigonometric operators in degrees.
> 
> This is what PostScript offers:
> 
> GS<2>dak@lola:/usr/local/tmp/lilypond/elisp$ gsnd
> GPL Ghostscript 9.15 (2014-09-22)
> Copyright (C) 2014 Artifex Software, Inc.  All rights reserved.
> This software comes with NO WARRANTY: see the file PUBLIC for details.
> GS>-1 0 atan ==
> 270.0
> GS>0 -1 atan ==
> 180.0
> GS>0 -1 atan sin ==
> 0.0
> GS>quit
> 
> This is also how METAFONT rolls:
> dak@lola:/usr/local/tmp/lilypond/elisp$ mf '\\relax'
> This is METAFONT, Version 2.7182818 (TeX Live 2015/dev/Debian) (preloaded
> base=mf)
> 
> *show angle(-1,0)
> 
> *;
> >> 180
> *show sin(angle(-1,0));
> >> sin
> ! Extra tokens will be flushed.
> <to be read again> 
>                    (
> <*> show sin(
>              angle(-1,0));
> ? 
> 
> *show sin;
> >> sin
> *show sind(angle(-1,0));
> >> 0
> *
> 
> If those graphical languages that we also work with decide it to be best to
keep
> off radians, maybe we should take a hint from them.

I can't follow all of your explanations, but it's clear what you demonstrate.

> 
> Our degree engine might even end up quite similar to what we have here, just
> working for multiples of 45 rather than PI/4, whether implemented in Scheme or
> C.  But the major difference is that the computer is _guaranteed_ to have an
> exact representation of multiples of 45.  In contrast, we have no guarantee
that
> 3*(PI/4) has an exact representation: this will only be the case if the lowest
> two significant bits of the representation of PI/4 are zero since 4 > 3pi/4 >
2
> and 1 > pi/4 > 0.5 and thus the mantissa of 3pi/4 takes up two more bits in
the
> most significant bits.
> 
> So we can't hope to keep stuff exact when juggling with multiples of angles in
> radians.
> 
> So let's rather discuss which degree-based primitives we'd rather have,
> implement them, and replace everything in LilyPond juggling with angles with
> them and make sure that our own commands only take degrees rather than
radians.
> 
> But I fear that this radian-based approach is only half a bandaid that will
> result in this problem reappearing in different forms.

Well, I'd happily agree, though I've not the slightest clue how to get (sin
what-ever-degree) work.
(You know I'm a schemer. If other languages come into play, I'm at a loss)

Please correct me, if I'm somewhere (or general) wrong:
All our trigonometric stuff (and guile itself) relies on radians, as an example:
`sin' is a guile-procedure _expecting_ a radians.
Thus, we would need to define our own sin-ly (or however we wold call it) and I
don't even know where to start such a task.

I thought about using polar-coordinates, but at some point PI (and not a degree)
needs to be introduced, see our `polar->rectangular'

Otoh, I'd really like a degree-based approach, any hint would be welcome.

[thomasmorley651, 2015-10-20 22:54:17]

no need for code-duplications, increasing readability of cyclic-base-value

[thomasmorley651, 2015-10-20 22:57:49]

On 2015/10/20 22:54:17, thomasmorley651 wrote:
> no need for code-duplications, increasing readability of cyclic-base-value

patch set 2 does some clean up, but not adressing comment #2

[thomasmorley651, 2015-10-20 23:00:47]

more clean up

[thomasmorley651, 2015-10-20 23:06:58]

https://codereview.appspot.com/269530043/diff/40001/scm/lily-library.scm
File scm/lily-library.scm (right):

https://codereview.appspot.com/269530043/diff/40001/scm/lily-library.scm#newc...
scm/lily-library.scm:751: (cyclic-base-value (+ value cycle) cycle))
what are those red dots??
Never seen them before.

[dak, 2015-10-21 09:32:07]

https://codereview.appspot.com/269530043/diff/40001/scm/lily-library.scm
File scm/lily-library.scm (right):

https://codereview.appspot.com/269530043/diff/40001/scm/lily-library.scm#newc...
scm/lily-library.scm:721: (moved-angle (- angle (/ (* quadrant PI) 2)))
Ok, let's just tackle this as a theoretical exercise since I remain unconvinced
that we are doing ourselves much of a favor straying from degrees except at the
final step.  Then angle-0-2pi is basically an independent source of numerical
ambiguity we want to avoid.  So I'd probably do it as follows:

(let* ((quadrant
	(inexact->exact (round (/ angle PI/2))))
       (moved-angle (- angle (* quadrant PI/2)))
       (s (sin moved-angle))
       (c (cos moved-angle))
       (x (coord-x coordinate))
       (y (coord-y coordinate)))
  (case (modulo quadrant 4)
    ((0) % -45 .. 45
     (cons ...))
    ((1) % 45 .. 135
     (cons ...))
    ((2) % 135 .. 225
     (cons ...))
    ((3) % 225 .. 315
     (cons ...))))

Note that this retains the problem that (* quadrant PI/2) might just
not be properly representable which can make the whole exercise end up
still not delivering perfect values for 3*PI/2 or 5*PI/2 or similar,
depending on just how those values were arrived at.  However,

(* exact-integer PI/2) likely has about the best chance to match
values that the user arrived with using similar expressions of his
own.  I'd also don't make two separate case statements for cos and sin
since that's likely a construct that is a bit expensive.

https://codereview.appspot.com/269530043/diff/40001/scm/lily-library.scm#newc...
scm/lily-library.scm:751: (cyclic-base-value (+ value cycle) cycle))
On 2015/10/20 23:06:58, thomasmorley651 wrote:
> what are those red dots??
> Never seen them before.

Spaces before tabs would be my guess.  Or just generally a spacing change. 
either of the two.

[thomasmorley651, 2015-10-25 13:23:38]

adressing David's comment

[thomasmorley651, 2015-10-25 13:27:38]

On 2015/10/21 09:32:07, dak wrote:
> https://codereview.appspot.com/269530043/diff/40001/scm/lily-library.scm
> File scm/lily-library.scm (right):
> 
>
https://codereview.appspot.com/269530043/diff/40001/scm/lily-library.scm#newc...
> scm/lily-library.scm:721: (moved-angle (- angle (/ (* quadrant PI) 2)))
> Ok, let's just tackle this as a theoretical exercise since I remain
unconvinced
> that we are doing ourselves much of a favor straying from degrees except at
the
> final step.  Then angle-0-2pi is basically an independent source of numerical
> ambiguity we want to avoid.  So I'd probably do it as follows:
> 
> (let* ((quadrant
> 	(inexact->exact (round (/ angle PI/2))))
>        (moved-angle (- angle (* quadrant PI/2)))
>        (s (sin moved-angle))
>        (c (cos moved-angle))
>        (x (coord-x coordinate))
>        (y (coord-y coordinate)))
>   (case (modulo quadrant 4)
>     ((0) % -45 .. 45
>      (cons ...))
>     ((1) % 45 .. 135
>      (cons ...))
>     ((2) % 135 .. 225
>      (cons ...))
>     ((3) % 225 .. 315
>      (cons ...))))
> 
> Note that this retains the problem that (* quadrant PI/2) might just
> not be properly representable which can make the whole exercise end up
> still not delivering perfect values for 3*PI/2 or 5*PI/2 or similar,
> depending on just how those values were arrived at.  However,
> 
> (* exact-integer PI/2) likely has about the best chance to match
> values that the user arrived with using similar expressions of his
> own.  I'd also don't make two separate case statements for cos and sin
> since that's likely a construct that is a bit expensive.

Changed it following your proposal.
Up to now I found no drawback.
> 
>
https://codereview.appspot.com/269530043/diff/40001/scm/lily-library.scm#newc...
> scm/lily-library.scm:751: (cyclic-base-value (+ value cycle) cycle))
> On 2015/10/20 23:06:58, thomasmorley651 wrote:
> > what are those red dots??
> > Never seen them before.
> 
> Spaces before tabs would be my guess.  Or just generally a spacing change. 
> either of the two.

Must be spacing change then

[benko.pal, 2015-10-30 13:20:35]

sorry to chime in that late, but:
am I right that the problem is that we get the rotation matrix
cos a   -sin a
sin a    cos a
inexact?  and if so, how the inexactness is present?  one of the diagonals is
exactly +/-1 while the other is not exactly 0?
in that case I'd suggest (would have suggested) to perform a paranoid
normalization of the (cos a, sin a) vector, which e.g. normalizes (1, 1e-9) into
(1, 0).

[dak, 2015-10-30 13:46:18]

On 2015/10/30 13:20:35, benko.pal wrote:
> sorry to chime in that late, but:
> am I right that the problem is that we get the rotation matrix
> cos a   -sin a
> sin a    cos a
> inexact?  and if so, how the inexactness is present?  one of the diagonals is
> exactly +/-1 while the other is not exactly 0?
> in that case I'd suggest (would have suggested) to perform a paranoid
> normalization of the (cos a, sin a) vector, which e.g. normalizes (1, 1e-9)
into
> (1, 0).

cos x = 1 - x^2/2 + ...

So basically this "paranoid normalization" would suppress all low values for
approximately |x| < sqrt (2 ulp) where "ulp" is the relative precision of our
floats.

So that looks like quite a noticeable cutoff relying on numerics.  The advantage
would be that the code is simpler, but once the code gets out of the
clamped-down region, the precision would not be better.

I really think we should move everything to degrees except for the final
calculations.

[benko.pal, 2015-10-30 15:07:42]

2015-10-30 14:46 GMT+01:00  <dak@gnu.org>:
> On 2015/10/30 13:20:35, benko.pal wrote:
>>
>> sorry to chime in that late, but:
>> am I right that the problem is that we get the rotation matrix
>> cos a   -sin a
>> sin a    cos a
>> inexact?  and if so, how the inexactness is present?  one of the
>
> diagonals is
>>
>> exactly +/-1 while the other is not exactly 0?
>> in that case I'd suggest (would have suggested) to perform a paranoid
>> normalization of the (cos a, sin a) vector, which e.g. normalizes
>> (1, 1e-9) into (1, 0).
>
>
> cos x = 1 - x^2/2 + ...
>
> So basically this "paranoid normalization" would suppress all low values
> for approximately |x| < sqrt (2 ulp) where "ulp" is the relative
> precision of our floats.
> So that looks like quite a noticeable cutoff relying on numerics.

I'm not sure I get that.  is this the description of the old problem
from a new viewpoint or an argument that if cos x equals 1 with
machine precision, we may still want a non-zero sin x?

> The
> advantage would be that the code is simpler, but once the code gets out
> of the clamped-down region, the precision would not be better.

sqrt (2 ulp) is quite a wide region; I'd think getting out of it means
we are in the wild anyway.  or that we are switching between single
and double (quadruple...) precision -- that must surely be avoided.

> I really think we should move everything to degrees except for the final
> calculations.

I'm not against that, but don't really know what are the final
calculations.  compute each (co)sine by checking first whether the
argument is an integer multiple of 90°?  or to review all our
calculations and reformulate so that each has at most one sin/cos/atan
computation?  does the problem come from functions in a "too" modular
hierarchy calling each other, converting between angles and sines back
and forth?  or would simply clamping all (co)sines with absolute value
below sqrt (2 ulp) to zero solve the practical problems at hand?

p

[dak, 2015-10-30 15:32:45]

On 2015/10/30 15:07:42, benko.pal wrote:
> 2015-10-30 14:46 GMT+01:00  <mailto:dak@gnu.org>:

> > I really think we should move everything to degrees except for the final
> > calculations.
> 
> I'm not against that, but don't really know what are the final
> calculations.

I'd argue for creating functions sind, cosd, angle in analogy to METAFONT's
functions.  And we'd use them exclusively so it might make some sense to define
C++ versions of them and export them in order to get consistent results.

They'd map everything to small degree angles before converting to/from radians.

> compute each (co)sine by checking first whether the
> argument is an integer multiple of 90°?

If you fold every operation to +/- 45 degrees before going to radians, you don't
need to check anything.  Multiples of 90 are perfectly representable in degrees.

> or to review all our
> calculations and reformulate so that each has at most one sin/cos/atan
> computation?  does the problem come from functions in a "too" modular
> hierarchy calling each other, converting between angles and sines back
> and forth?

No, it comes from pi/4 and its multiples not being representable in our floating
point arithmetic.

> or would simply clamping all (co)sines with absolute value
> below sqrt (2 ulp) to zero solve the practical problems at hand?

It would hide the problem.  But you cannot, for example, check angles for being
straight, even if their whole history was straightforward.  If we can't tell the
computer what we want it to do in its own terms, there are just no "trivial"
transformations or manipulations.

And PostScript (and PDF and METAFONT) _do_ represent angles in degrees.  So it's
sort of silly that we cannot get an angle of 180 degrees straight into
PostScript without change.

[wl_gnu.org, 2015-10-30 19:54:55]

> I'd argue for creating functions sind, cosd, angle in analogy to
> METAFONT's functions.  And we'd use them exclusively so it might
> make some sense to define C++ versions of them and export them in
> order to get consistent results.

Maybe there's a library out there which we could use, probably even
GNU stuff?


     Werner

[dak, 2015-10-30 20:03:23]

Werner LEMBERG <wl@gnu.org> writes:

>> I'd argue for creating functions sind, cosd, angle in analogy to
>> METAFONT's functions.  And we'd use them exclusively so it might
>> make some sense to define C++ versions of them and export them in
>> order to get consistent results.
>
> Maybe there's a library out there which we could use, probably even
> GNU stuff?

I'd be surprised as that is rather basic stuff.  Wait.

    pango_matrix_rotate ()

    void
    pango_matrix_rotate (PangoMatrix *matrix,
                         double degrees);

    Changes the transformation represented by matrix to be the
    transformation given by first rotating by degrees degrees
    counter-clockwise then applying the original transformation.

We are linking Pango anyway.  Maybe we should convert all our stencil
stuff into native Pango?  The unwholesome mess we pass around for path
expressions and convert back and forth is not particularly efficient.

And being all-Pango would make it easier to create more backends.  Even
if Pango's own PDF backend is not good enough (for dealing with
metadata, tables of content and similar stuff), the intermediate
representations should still be a good common denominator.

-- 
David Kastrup

[benko.pal, 2015-10-30 21:46:55]

> And PostScript (and PDF and METAFONT) _do_ represent angles in degrees.
> So it's sort of silly that we cannot get an angle of 180 degrees
> straight into PostScript without change.

now that convinced me.  if we want to output angles in degrees, we
don't want to switch back and forth neither between degrees and
radians, nor between angles and sines.  at most one atan (acos) call
somewhere.

p

[thomasmorley651, 2015-10-30 22:59:54]

On 2015/10/30 21:46:55, benko.pal wrote:
> > And PostScript (and PDF and METAFONT) _do_ represent angles in degrees.
> > So it's sort of silly that we cannot get an angle of 180 degrees
> > straight into PostScript without change.
> 
> now that convinced me.  if we want to output angles in degrees, we
> don't want to switch back and forth neither between degrees and
> radians, nor between angles and sines.  at most one atan (acos) call
> somewhere.
> 
> p

Just back home, noticing the discussion here.

Current patch tries to solve (or better: workaround) the problem in guile.
I have some other pending patches _needing_ exact values from coord-rotate.
I'd suggest to push it as is (unless some serious defect would be detected)
after countdown. So that I can go ahead with the follow up stuff.

Opinions?

Nevertheless, I'd love to work with degrees, rather than radians.
Though, I can't help implementing it as soon as C++ etc is involved.

[lemzwerg, 2015-10-31 05:27:58]

+1

[wl_gnu.org, 2015-10-31 05:33:44]

> Maybe we should convert all our stencil stuff into native Pango?

Sounds sensible.  However, I would add this to a TODO list than rather
working on it immediately, except if there are pressing needs.


    Werner

[dak, 2015-10-31 11:51:11]

On 2015/10/30 22:59:54, thomasmorley651 wrote:
> On 2015/10/30 21:46:55, benko.pal wrote:
> > > And PostScript (and PDF and METAFONT) _do_ represent angles in degrees.
> > > So it's sort of silly that we cannot get an angle of 180 degrees
> > > straight into PostScript without change.
> > 
> > now that convinced me.  if we want to output angles in degrees, we
> > don't want to switch back and forth neither between degrees and
> > radians, nor between angles and sines.  at most one atan (acos) call
> > somewhere.
> > 
> > p
> 
> Just back home, noticing the discussion here.
> 
> Current patch tries to solve (or better: workaround) the problem in guile.
> I have some other pending patches _needing_ exact values from coord-rotate.
> I'd suggest to push it as is (unless some serious defect would be detected)
> after countdown. So that I can go ahead with the follow up stuff.
> 
> Opinions?

I hope to get to a "degree rework", but this patch is sort of unrelated even
though I hope to eventually make it unneeded.  So it might as well go ahead for
now.

LGTM

[benko.pal, 2015-10-31 18:15:32]

> Just back home, noticing the discussion here.
>
> Current patch tries to solve (or better: workaround) the problem in
> guile.
> I have some other pending patches _needing_ exact values from
> coord-rotate.
> I'd suggest to push it as is (unless some serious defect would be
> detected) after countdown. So that I can go ahead with the follow up
> stuff.
>
> Opinions?

Fine with me.

> Nevertheless, I'd love to work with degrees, rather than radians.
> Though, I can't help implementing it as soon as C++ etc is involved.

I may help with c++.