@@ -301,49 +301,48 @@ jtQfromD(J jt, A w, void *yv, I mode) {
301301
302302static B
303303jtDfromQ (J jt, A w, void *yv) {
304- D d, f, n, *x, xb = (D)XBASE;
305- I cn, i, k, m, nn, pn, qn, r, *v, wn;
306- Q *wv;
307- X c, p, q, x2 = 0 ;
308- wn = AN (w);
309- wv = QAV (w);
310- x = (D *)yv;
311- nn = 308 / XBASEN;
312- for (i = 0 ; i < wn; ++i) {
313- p = wv[i].n ;
314- pn = AN (p);
315- k = 1 == pn ? AV (p)[0 ] : 0 ;
316- q = wv[i].d ;
317- qn = AN (q);
304+ auto const xb = (D)XBASE;
305+ auto const wn = AN (w);
306+ auto const wv = QAV (w);
307+ auto const x = (D *)yv;
308+ auto const nn = 308 / XBASEN;
309+
310+ // TODO: figure out nice algorithm for this
311+ auto const add_digits = [&](auto n, auto v) {
312+ auto f = 1.0 ;
313+ auto d = 0.0 ;
314+ DO (n, d += f * v[i]; f *= xb;);
315+ return d;
316+ };
317+
318+ X x2 = 0 ;
319+ for (int64_t i = 0 ; i < wn; ++i) {
320+ auto const p = wv[i].n ;
321+ auto const pn = AN (p);
322+ auto const k = 1 == pn ? AV (p)[0 ] : 0 ;
323+ auto const q = wv[i].d ;
324+ auto const qn = AN (q);
318325 if (k == XPINF)
319326 x[i] = inf;
320327 else if (k == XNINF)
321328 x[i] = infm;
322329 else if (pn <= nn && qn <= nn) {
323- n = 0.0 ;
324- f = 1.0 ;
325- v = AV (p);
326- DO (pn, n += f * v[i]; f *= xb;);
327- d = 0.0 ;
328- f = 1.0 ;
329- v = AV (q);
330- DO (qn, d += f * v[i]; f *= xb;);
331- x[i] = n / d;
330+ auto const n = add_digits (pn, AV (p));
331+ auto const d = add_digits (qn, AV (q));
332+ x[i] = n / d;
332333 } else {
333- k = 5 + qn;
334334 if (!x2)
335335 if (!(x2 = jtxc (jt, 2L ))) return 0 ;
336- if (!(c = jtxdiv (jt, jttake (jt, jtsc (jt, -(k + pn)), p), q, XMFLR))) return 0 ;
337- cn = AN (c);
338- m = MIN (cn, 5 );
339- r = cn - (m + k);
340- v = AV (c) + cn - m;
341- n = 0.0 ;
342- f = 1.0 ;
343- DO (m, n += f * v[i]; f *= xb;);
344- d = 1.0 ;
345- DQ (ABS (r), d *= xb;);
346- x[i] = 0 > r ? n / d : n * d;
336+ auto const k = 5 + qn;
337+ auto c = jtxdiv (jt, jttake (jt, jtsc (jt, -(k + pn)), p), q, XMFLR);
338+ if (!c) return 0 ;
339+ auto const cn = AN (c);
340+ auto const m = MIN (cn, 5 );
341+ auto const r = cn - (m + k);
342+ auto const v = AV (c) + cn - m;
343+ auto const n = add_digits (m, v);
344+ auto d = std::pow (xb, std::abs (r));
345+ x[i] = 0 > r ? n / d : n * d;
347346 }
348347 }
349348 return 1 ;
0 commit comments