Raul,
I thought it was just round off error. I noticed the same thing
when I rotate the vectors and then rotate them back:
rr=. vv qRot (uu;ww)
vv2=. rr qRot (uu;-ww)
the maximum difference, >./|vv2-vv, is about
~2e_15, while the smallest vv is ~1e_7. Most
of the differences are zero.
I've not worried about what the round off should
be. I've run programs with this method for days and
never noticed anything weird.
Cheers,
Patrick
On Wed, 13 Jul 2016, Raul Miller wrote:
I am a little concerned about the stability of this approach.
I extracted your timing test from earlier, and then started playing
with normalization.
rand=: [: ? ] # 0"_
vv=. 1e6 3$ _1+2*rand 3e6
uu=. norm 1e6 3$ _1+2*rand 3e6
ww=. _1p1+2p1*rand 1e6
echo timespacex'rr=. vv qRot (uu;ww)'
echo timespacex'rr2=. vv qRot ((norm uu);ww)'
assert rr-:rr2
This gives me:
0.232836 2.18107e8
0.274978 2.51662e8
|assertion failure: assert
| assert rr-:rr2
Subtracting I see a number of values with magnitude approximately
1.11022e_16, and I am not sure what to make of that.
What are the implications of norm being slightly unstable?
Thanks,
--
Raul
On Wed, Jul 13, 2016 at 12:42 AM, J. Patrick Harrington
<[email protected]> wrote:
That was even worse. All the LF's or CR's being removed.
I have been loading a file -- here's one last try, typed by hand :-(
-------------------------------------------------------------
qRot=: dyad define
'u w'=. y
u=. u*sin-:w
t=. +: u cross x
x+ (t*cos-:w)+ u cross t
)
L=: +/"1 &. (*:"_)
norm=: ] % L
cross=: (([: 1&|."1 [)* [: 2&|."1 ])- (([: 1&|."1 [)* [: 2&|."1 ])~
--------------------------------------------------------------------
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
