Hmm... I guess a part of my problem was figuring out how this was correct for very small values. The condition where it's executed is (pn<=nn&&qn<=nn)
Here, nn means that the term exceeds roughly 308 digits, and qn represents the denominator. And... testing shows that this actually doesn't work for very small values. A=: (?10x^350) % ?10x^50 B=: (?10x^50) % ?10x^350 0.0+A,B 2.72349e299 0 (0.0+(B*10x^50))%10^50 1.97574e_300 So.. I only need to focus on this for its significance when A is very large. Anyways, your perspective helped. A lot. Thanks, -- Raul On Thu, Oct 13, 2022 at 2:42 PM Henry Rich <henryhr...@gmail.com> wrote: > > 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 <rauldmil...@gmail.com> 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 ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
