160
:h2.DLdiv() - Divide Two Longs by One Long, Yielding Two Longs
162
Both the dividend and the divisor must be positive.
166
DLdiv(doublelong *quotient, /* also where dividend is, originally */
167
unsigned long divisor)
170
/* printf("DLdiv(%lx %lx)\n", quotient, divisor); */
171
*quotient /= divisor;
172
/* printf("DLdiv returns %lx\n", *quotient); */
174
register unsigned long u1u2 = quotient->high;
175
register unsigned long u3u4 = quotient->low;
176
register long u3; /* single digit of dividend */
177
register int v1,v2; /* divisor in registers */
178
register long t; /* signed copy of u1u2 */
179
register int qhat; /* guess at the quotient digit */
180
register unsigned long q3q4; /* low two digits of quotient */
181
register int shift; /* holds the shift value for normalizing */
182
register int j; /* loop variable */
184
/* printf("DLdiv(%x %x, %x)\n", quotient->high, quotient->low, divisor); */
186
* Knuth's algorithm works if the dividend is smaller than the
187
* divisor. We can get to that state quickly:
189
if (u1u2 >= divisor) {
190
quotient->high = u1u2 / divisor;
196
if (divisor <= MAXSHORT) {
199
* This is the case where the divisor is contained in one
200
* 'short'. It is worthwhile making this fast:
202
u1u2 = ASSEMBLE(u1u2, HIGHDIGIT(u3u4));
203
q3q4 = u1u2 / divisor;
205
u1u2 = ASSEMBLE(u1u2, LOWDIGIT(u3u4));
206
quotient->low = ASSEMBLE(q3q4, u1u2 / divisor);
212
* At this point the divisor is a true 'long' so we must use
215
* Step D1: Normalize divisor and dividend (this makes our 'qhat'
216
* guesses more accurate):
218
for (shift=0; !SIGNBITON(divisor); shift++, divisor <<= 1) { ; }
222
if ((u1u2 >> (LONGSIZE - shift)) != 0 && shift != 0)
223
Abort("DLdiv: dividend too large");
224
u1u2 = (u1u2 << shift) + ((shift == 0) ? 0 : u3u4 >> (LONGSIZE - shift));
228
* Step D2: Begin Loop through digits, dividing u1,u2,u3 by v1,v2,
229
* then shifting U left by 1 digit:
231
v1 = HIGHDIGIT(divisor);
232
v2 = LOWDIGIT(divisor);
234
u3 = HIGHDIGIT(u3u4);
236
for (j=0; j < 2; j++) {
239
* Step D3: make a guess (qhat) at the next quotient denominator:
241
qhat = (HIGHDIGIT(u1u2) == v1) ? MAXSHORT : u1u2 / v1;
243
* At this point Knuth would have us further refine our
244
* guess, since we know qhat is too big if
246
* v2 * qhat > ASSEMBLE(u1u2 % v, u3)
248
* That would make sense if u1u2 % v was easy to find, as it
249
* would be in assembly language. I ignore this step, and
250
* repeat step D6 if qhat is too big.
254
* Step D4: Multiply v1,v2 times qhat and subtract it from
259
* The high digit of u3 now contains the "borrow" for the
260
* rest of the substraction from u1,u2.
261
* Sometimes we can lose the sign bit with the above.
262
* If so, we have to force the high digit negative:
266
t |= -1 << SHORTSIZE;
267
t += u1u2 - qhat * v1;
268
/* printf("..>divide step qhat=%x t=%x u3=%x u1u2=%x v1=%x v2=%x\n",
269
qhat, t, u3, u1u2, v1, v2); */
270
while (t < 0) { /* Test is Step D5. */
273
* D6: Oops, qhat was too big. Add back in v1,v2 and
274
* decrease qhat by 1:
276
u3 = LOWDIGIT(u3) + v2;
277
t += HIGHDIGIT(u3) + v1;
279
/* printf("..>>qhat correction t=%x u3=%x qhat=%x\n", t, u3, qhat); */
282
* Step D7: shift U left one digit and loop:
285
if (HIGHDIGIT(u1u2) != 0)
286
Abort("divide algorithm error");
287
u1u2 = ASSEMBLE(u1u2, LOWDIGIT(u3));
289
q3q4 = ASSEMBLE(q3q4, qhat);
291
quotient->low = q3q4;
292
/* printf("DLdiv returns %x %x\n", quotient->high, quotient->low); */
298
:h3.DLadd() - Add Two Double Longs
300
In this case, the doublelongs may be signed. The algorithm takes the
301
piecewise sum of the high and low longs, with the possibility that the
302
high should be incremented if there is a carry out of the low. How to
303
tell if there is a carry? Alex Harbury suggested that if the sum of
304
the lows is less than the max of the lows, there must have been a
305
carry. Conversely, if there was a carry, the sum of the lows must be
306
less than the max of the lows. So, the test is "if and only if".
310
DLadd(doublelong *u, /* u = u + v */
314
/* printf("DLadd(%lx %lx)\n", *u, *v); */
316
/* printf("DLadd returns %lx\n", *u); */
318
register unsigned long lowmax = MAX(u->low, v->low);
320
/* printf("DLadd(%x %x, %x %x)\n", u->high, u->low, v->high, v->low); */
328
:h3.DLsub() - Subtract Two Double Longs
330
Testing for a borrow is even easier. If the v.low is greater than
331
u.low, there must be a borrow.
335
DLsub(doublelong *u, /* u = u - v */
339
/* printf("DLsub(%lx %lx)\n", *u, *v); */
341
/* printf("DLsub returns %lx\n", *u); */
343
/* printf("DLsub(%x %x, %x %x)\n", u->high, u->low, v->high, v->low);*/
351
160
:h3.DLrightshift() - Macro to Shift Double Long Right by N
401
210
return ((negative) ? -ret : ret);
406
:h3.FPdiv() - Divide Two Fractional Pel Values
408
These values may be signed. The function returns the quotient.
412
FPdiv(fractpel dividend, fractpel divisor)
414
doublelong w; /* result will be built here */
415
int negative = FALSE; /* flag for sign bit */
417
register fractpel ret;
421
dividend = -dividend;
426
negative = !negative;
429
w.low = dividend << FRACTBITS;
430
w.high = dividend >> (LONGSIZE - FRACTBITS);
432
if (w.high != 0 || SIGNBITON(w.low)) {
433
w.low = TOFRACTPEL(MAXSHORT);
435
return( (negative) ? -w.low : w.low);
437
w = ((long)dividend) << FRACTBITS;
439
if (w & 0xffffffff80000000L ) {
440
ret = TOFRACTPEL(MAXSHORT);
444
return( (negative) ? -ret : ret);
449
:h3.FPstarslash() - Multiply then Divide
451
Borrowing a chapter from the language Forth, it is useful to define
452
an operator that first multiplies by one constant then divides by
453
another, keeping the intermediate result in extended precision.
457
FPstarslash(fractpel a, /* result = a * b / c */
461
doublelong w; /* result will be built here */
462
int negative = FALSE;
464
register fractpel ret;
467
if (a < 0) { a = -a; negative = TRUE; }
468
if (b < 0) { b = -b; negative = !negative; }
469
if (c < 0) { c = -c; negative = !negative; }
474
if (w.high != 0 || SIGNBITON(w.low)) {
475
w.low = TOFRACTPEL(MAXSHORT);
477
return((negative) ? -w.low : w.low);
479
if (w & 0xffffffff80000000L ) {
480
ret = TOFRACTPEL(MAXSHORT);
484
return( (negative) ? -ret : ret);
added
removed
src/Type1/arith.c
