kicad-developers team mailing list archive

[Lorenzo Marcantonio <l.marcantonio@xxxxxxxxxxxx>]

We talked about this the other day (when discussing the fully fixed
point implementation of the CERN geometric lib). I think that many
things with angle can't be 'detached' from floating point.

Rotation by some angle is one of these, sin and cos are ugly beasts
(from a mathematical stand point); well, whereever appears PI or sqrt(2)
or (save us from it:D) e, there *will* be some numerical trouble.
I think that even using some fixed point implementation of sin and cos
there wouldn't be much to be gained on the stability side (correct me if
I'm wrong, I took numerical methods too many years ago).

So when I say: rotate this pad 45 degrees, there is no hope. However 45
is constant so the sin(45) and cos(45) values are stable, wherever you
put them (in fact, in ancient kicad, they where tabulated!)

However a common pattern I've seen more than twice is the following:
- Take a vector (or two points)
- Compute the angle between them
- Use this angle to rotate other stuff

Here we have: x,y (as int) 
-> x,y (as double, 'theorically' no loss at all) 
-> y/x (things begin to get ugly if you know how a divide works...)
-> alpha = atan(y/x) (power series, most probably)
-> (Usually) alpha multiplied to go from radian to degrees
-> (Usually) alpha multiplied to go back
-> sin(alpha), cos(alpha) (another two power series, in fact the x87 can
                           compute both at the same time, maybe gcc is
                           smart enough to notice it since these are
                           intrinsics :D)
-> Transform the point using them (let's say four multiply and four add)

If we're lucky the compiler will leave the stuff on the FPU stack so we
get long double precision. Also alpha would be often some nasty number
like PI/4, however seems that the FPU (or the math libs) handle this
precisely enough. Also the atan2 could be so smart to check for ugly
divisions and do acot(x/y) instead of atan(y/x). Or maybe the double
precision is simply enough for all our needs (16 digit against 9 and
a bit of the 32 bit integer).

However (this is basic linear algebra) we can see that:
- Taken the original vector, we normalize it.
- By definition, the resulting unity vector is the cos(alpha),sin(alpha)
  pair we want
- Transform the point like before

Looking at the new series of operations:
-> x,y (to double exactly as before)
-> need the norm: two squares, a sum and a square root
   (I have no idea if the sqrt is done with a series or some dedicated
   algorithm, however I 'feel' that's 'easier' on the number than an atan)
-> Two division to get cos(alpha) and sin(alpha) (here the division
   could go nasty only if the vector is *extremely* short; since we're
   starting with integer coordinated I don't think that could happen;
   correct me if I'm wrong)
-> Point transform exactly like before

In this way we could avoid a (risky) division and *three* trigonometric
calls, in exchange for a square root and a two (relatively safe)
division.

However this is mostly gut feeling, I've no real experience about
numerics with trascendents... can someone tell me if there would be
a *useful* gain in numeric stability in doing this? the idea would be
a RotatePoint family taking the reference vector instead of the angle in
decidegrees, obviously.

Since this is a somewhat 'strange' thing to do (even if actually easy to
do) I feel is not warranted if it has no numeric advantage. I already
did this in my push_track extension (extracting a sin from a dot
product) and the numbers where 'rounder' but probably this was only
a chance.

PS: Did this experiment in octave (matlab if you like). It's long but the
conclusions substantially are: until someone give me a better reason, there is
no need to do this; it seems that double trig is precise enough for what we are
doing. Had to try before saying this:D

format long;
octave:1> format long;
octave:2> pi/4
ans =  0.785398163397448
octave:3> atan2(10000, 10000)
ans =  0.785398163397448
octave:4> ans * 1800.0 / pi
ans =  450

OK, so far, so good. The sacred 45 degrees angle is handled correctly

octave:5> sin(pi/4)
ans =  0.707106781186547
octave:6> cos(pi/4)
ans =  0.707106781186548

(yes, I am cheating, but the last number is not the same:P:P:P)

octave:7> rot=[ cos(pi/4) -sin(pi/4); sin(pi/4) cos(pi/4)]
rot =

   0.707106781186548  -0.707106781186547
   0.707106781186547   0.707106781186548

octave:8> [ 12345678 23456789 ] * rot
ans =

   25316167.19890759    7856741.93461644

(it's not a serious test, it's just to get an idea;
That would be 25316167, 7856742 once rounded)

OK, now for method 2:

octave:10> sqrt(10000*10000+10000*10000)
ans =  14142.1356237310
octave:11> 10000/ans
ans =  0.707106781186548
(the theory is right, at least:D)

octave:12> rot = [ 0.707106781186548 -0.707106781186548;
0.707106781186548 0.707106781186548]
rot =

   0.707106781186548  -0.707106781186548
   0.707106781186548   0.707106781186548
octave:13> [ 12345678 23456789 ] * rot
ans =

   25316167.19890761    7856741.93461645

(yep, thats the same, for all the practical reasons. We're talking about
femtometers errors, here...)

So, for the 'main case', no advantage. I actually *hoped* it!
Now I'm curious of seeing if tan/atan is sufficiently smart (I think it
should, we only use half of the double significant digits)

octave:18> atan(1234567890)
ans =  1.57079632598490
octave:19> tan(1.57079632598490)
ans =  1234573120.67690

I tought better... an error on the 10µm position is still acceptable
(also it's the angle we need, not the tan itself)

octave:2> 1/1234567890
ans =  8.10000007371000e-10
octave:3> atan(ans)
ans =  8.10000007371000e-10
octave:4> tan(ans)
ans =  8.10000007371000e-10
octave:5> 

As expected, just like sin(x)=x for sufficiently small x, it survived
the trip.

octave:2> atan(1234567890/1234)
ans =  1.57079532725489
octave:3> acot(1234/1234567890)
ans =  1.57079532725489
octave:4> atan2(1234567890,1234)                                                            
ans =  1.57079532725489
octave:5> atan(1234/1234567890)
ans =  9.99540009095481e-07
octave:6> acot(1234567890/1234)
ans =  9.99540009095481e-07
octave:7> atan2(1234,1234567890)
ans =  9.99540009095481e-07

As expected there are sufficient working digits to do the work. At least
the angle come out the same however you compute it (obviously the
IEEE754 people had done they homework!)

Now, the important thing: does the cos/sin of the angle come out
correct? I don't care about tan being reversible, I want to arctan and
then sin the angle...

Note: the test vector are purposely orthogonal (even if a bit silly). As a
magnitude, they could happen, for example, intersecting nearly parallel lines
or simply with boards large than one meter:D also that's meant as a (trivial)
stress test (BTW using 1,1234567890 for atan2 pulls the cosine to 1, which I hope is acceptable)

octave:2> alpha =atan2(1234,1234567890)
alpha =  9.99540009095481e-07
octave:3> rotalpha=[cos(alpha) -sin(alpha);sin(alpha) -cos(alpha)]
rotalpha =

   9.99999999999501e-01  -9.99540009095315e-07
   9.99540009095315e-07  -9.99999999999501e-01

octave:4> beta=atan2(1234567890,1234)
beta =  1.57079532725489
octave:5> rotbeta=[cos(beta) -sin(beta);sin(beta) -cos(beta)]
rotbeta =

   9.99540009044568e-07  -9.99999999999501e-01
   9.99999999999501e-01  -9.99540009044568e-07

OK, the difference is on the 1e-17 point (if I've counted right), so it seems
there are no problems here. Now for the alternative method:

octave:2> nrm = sqrt(1234567890*1234567890 + 1234*1234)
nrm =  1234567890.00062
octave:3> [1234567890/nrm 1234/nrm]
ans =

   9.99999999999501e-01   9.99540009095315e-07

Still the error in the 1e-17. Does it *really* matter? I hope not, since it's a
relative error...

Last thing: this is a (relatively) huge nearly 45 degrees vector, at the limit of the int range; the squaring in fact would need more digits than the available ones)

octave:7> alpha =atan2(2134567890,2134565432)
alpha =  0.785398739158338
octave:8> rotalpha=[cos(alpha) -sin(alpha);sin(alpha) -cos(alpha)]
rotalpha =

   0.707106374062001  -0.707107188310860
   0.707107188310860  -0.707106374062001
octave:9> nrm=sqrt(2134567890*2134567890+2134565432*2134565432)                             
nrm =  3018733121.77615
octave:10> [2134567890/nrm 2134565432/nrm]
ans =

   0.707107188310860   0.707106374062001

It came out *exactly* the same!

I agree that this is *not* scientific proof, but at least now I'm confident that (in the int range!) double trig functions are pretty safe to use, at least in general. On Intel CPUs, of course (AFAIK Java reimplemented them because they where not agreeing with their specifications).

So, while interesting, I feel that (unless I forgot something) there is no need to change the way to do rotations.

-- 
Lorenzo Marcantonio
Logos Srl