J does not use long double.
It is possible that 8.07 used 8087 instructions with 80-bit
floating-point; I don't remember.
The QR algorithm in JE is a clever recursive algorithm that Roger
invented. It uses +/ . * to do almost all the calculations; reciprocals
are performed only at the bottom of the recursion, on 2x2 or 1x1 matrices.
At the lower levels of recursion it works with short, wide matrices that
the blocking matrix-multiply code can't grip efficiently, so I added a
version of +/ . * that does matrix multiplication via inner products two
at a time: just the thing for such cases. The overall algorithm is
unchanged.
The differences you report look pretty small to me, and could be
artifacts of differing orders of summation.
I have seen references to what Python does, but I don't know what its
qr() code is defined to do. If the matrix is ill-conditioned, perhaps
some extra processing is performed.
Henry Rich
On 5/2/2021 12:54 PM, Raul Miller wrote:
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,
--
This email has been checked for viruses by AVG.
https://www.avg.com
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
