The code you post divides (n*1e20)/d and treats the low 5 digits of result as a fractional part. It multiplies the fraction by 1e20 and then divides that by 1e20, rounding. (It may use fewer digits in some cases.). The rounded result is transferred to the quotient.
I don't see that this is guaranteed to be correct for big values. If you need to round the quotient, why don't you add half the denominator to the numerator before using the gmp function? You might need to multiply both parts by 2 first. Henry Rich On Thu, Oct 13, 2022, 10:49 AM Raul Miller <[email protected]> wrote: > 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 > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
