https://github.com/jsoftware/jsource/blob/master/jsrc/k.c#L306
I've been spinning my wheels on this, in the context of https://github.com/jsoftware/jsource/blob/master/test/gxco.ijs#L92 In other words, where X p is the numerator of a rational and I pn is the number of "big digits" in p and X q is the denominator of a rational and I qn is the number of "big digits" in q, and p and q are near the limits of values which can represented as IEEE 754 floating points (type D), k=5+qn; if(!(c=xdiv(take(sc(-(k+pn)),p),q,XMFLR)))R 0; cn=AN(c); m=MIN(cn,5); r=cn-(m+k); v=AV(c)+cn-m; n=0.0; f=1.0; DO(m, n+=f*v[i]; f*=xb;); d=1.0; DQ(ABS(r), d*=xb;); x[i]=0>r?n/d:n*d; somehow is more stable than any rewrite than I can come up with. (I keep getting off-by-one least significant bit errors in my rewrites). I think that part of this is that I do not adequately understand what xdiv(take(sc(-(k+pn)),p),q,XMFLR) accomplishes. Part of this, of course, is that I'm trying to come up with something which works where xb is something other than 1e4 (which would be very slow, for me). Anyway... here's an example of a J model of a bad rewrite: base=: x:1e4 base=: 2^64x basedigits=: >.10^.base X2D2=: {{ if. 0=y do. 0.0 else. yn=. >.base^.y yl=. 0.0+(yn#base)#:|y ys=. *y ys*base#.yl end. }} R2Db=: {{ y=. y%+./y assert. (,2) -: $y 'N D'=. y 'en ed'=. (>.base^.|y)-<.308%basedigits e=. 0 if. 0<en do. N=. N <.@% base^en [ e=. en end. if. 0<ed do. D=. D <.@% base^ed [ e=. e-ed end. (N %&X2D2 D)*<.base^e }} R2D=: R2Db"1 f1=: 4 : 0 p=: (_1^x ?@$ 2) * x ?@$ y q=: 2+($p) ?@$ 20 e=: (_1^($p) ?@$ 2),.q^x:<.350*q^.10 e=: e % +./"1|e d=: p - R2D p (({.@]+*&{:),.{:@])"0 1 e assert. 0 = d 1 ) 1000 f1"0 ] 10^2 3 This seems like it should be simple, but I've been staring at it too long and am not seeing the forest because of all the trees... If anyone has insight on this one, I'd love to hear it. Thanks, -- Raul ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
