arith.c

远程桌面连接工具

               */
 
               /*
               * Step D4: Multiply v1,v2 times qhat and subtract it from
               * u1,u2,u3:
               */
               u3 -= qhat * v2;
               /*
               * The high digit of u3 now contains the "borrow" for the
               * rest of the substraction from u1,u2.
               * Sometimes we can lose the sign bit with the above.
               * If so, we have to force the high digit negative:
               */
               t = HIGHDIGIT(u3);
               if (t > 0)
                       t |= -1 << SHORTSIZE;
               t += u1u2 - qhat * v1;
/* printf("..>divide step qhat=%x t=%x u3=%x u1u2=%x v1=%x v2=%x\n",
                             qhat, t, u3, u1u2, v1, v2); */
               while (t < 0) {  /* Test is Step D5.                          */
 
                       /*
                       * D6: Oops, qhat was too big.  Add back in v1,v2 and
                       * decrease qhat by 1:
                       */
                       u3 = LOWDIGIT(u3) + v2;
                       t += HIGHDIGIT(u3) + v1;
                       qhat--;
/* printf("..>>qhat correction t=%x u3=%x qhat=%x\n", t, u3, qhat); */
               }
               /*
               * Step D7:  shift U left one digit and loop:
               */
               u1u2 = t;
               if (HIGHDIGIT(u1u2) != 0)
                       abort("divide algorithm error");
               u1u2 = ASSEMBLE(u1u2, LOWDIGIT(u3));
               u3 = LOWDIGIT(u3u4);
               q3q4 = ASSEMBLE(q3q4, qhat);
       }
       quotient->low = q3q4;
/* printf("DLdiv returns %x %x\n", quotient->high, quotient->low); */
#endif /* !LONG64 */
       return;
}
 
/*
:h3.DLadd() - Add Two Double Longs
 
In this case, the doublelongs may be signed.  The algorithm takes the
piecewise sum of the high and low longs, with the possibility that the
high should be incremented if there is a carry out of the low.  How to
tell if there is a carry?  Alex Harbury suggested that if the sum of
the lows is less than the max of the lows, there must have been a
carry.  Conversely, if there was a carry, the sum of the lows must be
less than the max of the lows.  So, the test is "if and only if".
*/
 
void DLadd(u, v)
       doublelong *u;        /* u = u + v                                    */
       doublelong *v;
{
#ifdef LONG64
/* printf("DLadd(%lx %lx)\n", *u, *v); */
       *u = *u + *v;
/* printf("DLadd returns %lx\n", *u); */
#else
       register unsigned long lowmax = MAX(u->low, v->low);
 
/* printf("DLadd(%x %x, %x %x)\n", u->high, u->low, v->high, v->low); */
       u->high += v->high;
       u->low += v->low;
       if (lowmax > u->low)
               u->high++;
#endif
}
/*
:h3.DLsub() - Subtract Two Double Longs
 
Testing for a borrow is even easier.  If the v.low is greater than
u.low, there must be a borrow.
*/
 
void DLsub(u, v)
       doublelong *u;        /* u = u - v                                    */
       doublelong *v;
{
#ifdef LONG64
/* printf("DLsub(%lx %lx)\n", *u, *v); */
       *u = *u - *v;
/* printf("DLsub returns %lx\n", *u); */
#else
/* printf("DLsub(%x %x, %x %x)\n", u->high, u->low, v->high, v->low);*/
       u->high -= v->high;
       if (v->low > u->low)
               u->high--;
       u->low -= v->low;
#endif
}
/*
:h3.DLrightshift() - Macro to Shift Double Long Right by N
*/
 
/*SHARED LINE(S) ORIGINATED HERE*/
 
/*
:h2.Fractional Pel Arithmetic
*/
/*
:h3.FPmult() - Multiply Two Fractional Pel Values
 
This funtion first calculates w = u * v to "doublelong" precision.
It then shifts w right by FRACTBITS bits, and checks that no
overflow will occur when the resulting value is passed back as
a fractpel.
*/
 
fractpel FPmult(u, v)
  register fractpel u,v;
{
  doublelong w;
  register int negative = FALSE; /* sign flag */
#ifdef LONG64
  register fractpel ret;
#endif
 
  if ((u == 0) || (v == 0)) return (0);
 
 
  if (u < 0) {u = -u; negative = TRUE;}
  if (v < 0) {v = -v; negative = !negative;}
 
  if (u == TOFRACTPEL(1)) return ((negative) ? -v : v);
  if (v == TOFRACTPEL(1)) return ((negative) ? -u : u);
 
  DLmult(&w, u, v);
  DLrightshift(w, FRACTBITS);
#ifndef LONG64
  if (w.high != 0 || SIGNBITON(w.low)) {
        IfTrace2(TRUE,"FPmult: overflow, %px%p\n", u, v);
        w.low = TOFRACTPEL(MAXSHORT);
  }
 
  return ((negative) ? -w.low : w.low);
#else
  if (w & 0xffffffff80000000L ) {
        IfTrace2(TRUE,"FPmult: overflow, %px%p\n", u, v);
        ret = TOFRACTPEL(MAXSHORT);
  }
  else
        ret = (fractpel)w;
 
  return ((negative) ? -ret : ret);
#endif
}
 
/*
:h3.FPdiv() - Divide Two Fractional Pel Values
 
These values may be signed.  The function returns the quotient.
*/
 
fractpel FPdiv(dividend, divisor)
       register fractpel dividend;
       register fractpel divisor;
{
       doublelong w;         /* result will be built here                    */
       int negative = FALSE; /* flag for sign bit                            */
#ifdef LONG64
       register fractpel ret;
#endif
 
       if (dividend < 0) {
               dividend = -dividend;
               negative = TRUE;
       }
       if (divisor < 0) {
               divisor = -divisor;
               negative = !negative;
       }
#ifndef LONG64
       w.low = dividend << FRACTBITS;
       w.high = dividend >> (LONGSIZE - FRACTBITS);
       DLdiv(&w, divisor);
       if (w.high != 0 || SIGNBITON(w.low)) {
               IfTrace2(TRUE,"FPdiv: overflow, %p/%p\n", dividend, divisor);
               w.low = TOFRACTPEL(MAXSHORT);
       }
       return( (negative) ? -w.low : w.low);
#else
       w = ((long)dividend) << FRACTBITS;
       DLdiv(&w, divisor);
       if (w & 0xffffffff80000000L ) {
               IfTrace2(TRUE,"FPdiv: overflow, %p/%p\n", dividend, divisor);
               ret = TOFRACTPEL(MAXSHORT);
       }
       else
               ret = (fractpel)w;
       return( (negative) ? -ret : ret);
#endif
}
 
/*
:h3.FPstarslash() - Multiply then Divide
 
Borrowing a chapter from the language Forth, it is useful to define
an operator that first multiplies by one constant then divides by
another, keeping the intermediate result in extended precision.
*/
 
fractpel FPstarslash(a, b, c)
       register fractpel a,b,c;  /* result = a * b / c                       */
{
       doublelong w;         /* result will be built here                    */
       int negative = FALSE;
#ifdef LONG64
       register fractpel ret;
#endif
 
       if (a < 0) { a = -a; negative = TRUE; }
       if (b < 0) { b = -b; negative = !negative; }
       if (c < 0) { c = -c; negative = !negative; }
 
       DLmult(&w, a, b);
       DLdiv(&w, c);
#ifndef LONG64
       if (w.high != 0 || SIGNBITON(w.low)) {
               IfTrace3(TRUE,"FPstarslash: overflow, %p*%p/%p\n", a, b, c);
               w.low = TOFRACTPEL(MAXSHORT);
       }
       return((negative) ? -w.low : w.low);
#else
       if (w & 0xffffffff80000000L ) {
               IfTrace3(TRUE,"FPstarslash: overflow, %p*%p/%p\n", a, b, c);
               ret = TOFRACTPEL(MAXSHORT);
       }
       else
               ret = (fractpel)w;
       return( (negative) ? -ret : ret);
#endif
}