I think you should simply use the result of mpq_get_d. It truncates toward 0 which is ok i think.
I understand the code fragment a bit better now. It does the divide to get the integer part plus 60 or so fraction bits which it then adds to the integer part. You could get the same effect by adding 1 limb (instead of 5) and multiplying the extra value by the right power of 2. Henry Rich On Thu, Oct 13, 2022, 4:04 PM Raul Miller <rauldmil...@gmail.com> wrote: > 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 > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
