I don't understand the question: what are you trying to compute? Doesn't gmp support division?
Knuth has code for division IIRC. hhr 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
