Here's something interesting.
When I compare the 128!:0 implementation of QR with a straightforward
J implementation of QR, I get slightly different results.
QR0=: 128!:0
QR1=: 3 : 0
n=.{:$A=.y
if. 1>:n do.
A ((% {.@,) ; ]) %:(+@|: A) +/ .* A
else.
m =.>.n%2
A0=.m{."1 A
A1=.m}."1 A
'Q0 R0'=.QR1 A0
'Q1 R1'=.QR1 A1 - Q0 +/ .* T=.(+@|: Q0) +/ .* A1
(Q0,.Q1);(R0,.T),(-n){."1 R1
end.
)
In various J90x version (I tested J901, J902 and J903 beta):
(QR0 -&> QR1) p:i.3 3
0 _2.22045e_16 2.33147e_15
5.55112e_17 _8.88178e_16 9.74221e_15
1.11022e_16 _2.16493e_15 _6.96318e_15
0 3.55271e_15 7.10543e_15
0 1.33227e_15 _1.5099e_14
0 0 _2.22045e_16
x=:_0.5+(i.1000)%999
c=:(^x) %. X=:(^ 0j1 * x */ i:10)
(+/%#)*:|(X (+/ .*) c)-(^x)
6.57613e_5
Compare J807:
(QR0 -&> QR1) p:i.3 3
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
x=:_0.5+(i.1000)%999
c=:(^x) %. X=:(^ 0j1 * x */ i:10)
(+/%#)*:|(X (+/ .*) c)-(^x)
0.0279006
I imagine that this has to do with floating point precision
limitations, though I have not looked deeply enough to characterize
the issues here. That said, when I look at the current J source, I see
that it is using 'long double' (which has a 64 bit mantissa):
jsource$ find . -iname '*.h' | xargs grep -l 'long double' | xargs
grep -c 'long double'
./jsrc/jtype.h:1
./sleef/include/sleef.h:4
./sleef/src/arch/helpervecext.h:15
./sleef/src/common/misc.h:4
So that's probably part of it.
FYI,
--
Raul
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
