Matrix Division on the 28C/S (and later)

[Werner]

Hello Thomas, no that is correct, the FMA 'compensated dot' algorithm is actually not the same as just using 68 digits ;-) As long as we're just performing a DOT product with one argument being a vector with all 1's, we're performing additions only, and we can actually do the same in RPN:
(based on Dekker's double precision algorithms, but modified a bit to work for decimal floating point arithmetic):

Code:

extended precision summation
In: [x]
Out: Sum(x[i])

>LBL "\S"
 EDIT
 ABS
 CLST
 LASTX
>LBL 10
 XEQ DS+
 J+
 RCLEL
 FC? 77
 GTO 10
 EXITALL
 CLX
 +
 + @ leave this off to see the full double-precision result
 END

@ add a number to a double-precision number
@        L    X    Y    Z    T
>LBL "DS+"        x    ts    hs
 X<>Y
 X<> ST Z        hs    x    ts
 XEQ "S+"        tr    hr    ts
 X<>Y
 X<> ST Z        ts    tr    hr
 +
>LBL "Q+" @ 'quick add 2 numbers'
 X<>Y
 RCL+ ST Y
 STO- ST L
 X<>Y
 RCL+ ST L
 RTN
>LBL "S+" @ add 2 numbers, result in double precision
 LSTO "."
 RCL+ ST Y
 RCL- ST Y
 STO- "."
 X<> ST L
 STO- ST L
 X<>Y
 RCL+ ST L
 RCL+ "."
 END

I can probably do the same for the full-fledged DOT.. I have the 'multiply two numbers and obtain a double-precision result' routines as well. The only thing I don't have for now is time ;-)

Cheers, Werner