Back to index

tetex-bin  3.0
Defines | Functions
arith.c File Reference
#include <stdio.h>
#include "types.h"
#include "objects.h"
#include "spaces.h"
#include "arith.h"

Go to the source code of this file.

Defines

#define ASSEMBLE(hi, lo)   ((((ULONG)hi)<<SHORTSIZE)+(lo))
#define HIGHDIGIT(u)   ((u)>>SHORTSIZE)
#define LOWDIGIT(u)   ((u)&MAXSHORT)
#define SIGNBITON(w)   (((LONG)w)<0)

Functions

void DLmult (doublelong *product, ULONG u, ULONG v)
void DLdiv (doublelong *quotient, ULONG divisor)
void DLadd (doublelong *u, doublelong *v)
void DLsub (doublelong *u, doublelong *v)
fractpel FPmult (fractpel u, fractpel v)
fractpel FPdiv (fractpel dividend, fractpel divisor)
fractpel FPstarslash (fractpel a, fractpel b, fractpel c)

Define Documentation

#define ASSEMBLE (   hi,
  lo 
)    ((((ULONG)hi)<<SHORTSIZE)+(lo))

Definition at line 87 of file arith.c.

#define HIGHDIGIT (   u)    ((u)>>SHORTSIZE)

Definition at line 92 of file arith.c.

#define LOWDIGIT (   u)    ((u)&MAXSHORT)

Definition at line 93 of file arith.c.

#define SIGNBITON (   w)    (((LONG)w)<0)

Definition at line 98 of file arith.c.


Function Documentation

void DLadd ( doublelong u,
doublelong v 
)

Definition at line 293 of file arith.c.

{
       register ULONG lowmax = TYPE1_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++;
}
void DLdiv ( doublelong quotient,
ULONG  divisor 
)

Definition at line 155 of file arith.c.

{
       register ULONG u1u2 = quotient->high;
       register ULONG u3u4 = quotient->low;
       register LONG u3;     /* single digit of dividend                     */
       register int v1,v2;   /* divisor in registers                         */
       register LONG t;      /* signed copy of u1u2                          */
       register int qhat;    /* guess at the quotient digit                  */
       register ULONG q3q4;  /* low two digits of quotient           */
       register int shift;   /* holds the shift value for normalizing        */
       register int j;       /* loop variable                                */
 
/* printf("DLdiv(%x %x, %x)\n", quotient->high, quotient->low, divisor); */
       /*
       * Knuth's algorithm works if the dividend is smaller than the
       * divisor.  We can get to that state quickly:
       */
       if (u1u2 >= divisor) {
               quotient->high = u1u2 / divisor;
               u1u2 %= divisor;
       }
       else
               quotient->high = 0;
 
       if (divisor <= MAXSHORT) {
 
               /*
               * This is the case where the divisor is contained in one
               * 'short'.  It is worthwhile making this fast:
               */
               u1u2 = ASSEMBLE(u1u2, HIGHDIGIT(u3u4));
               q3q4 = u1u2 / divisor;
               u1u2 %= divisor;
               u1u2 = ASSEMBLE(u1u2, LOWDIGIT(u3u4));
               quotient->low = ASSEMBLE(q3q4, u1u2 / divisor);
               return;
       }
 
 
       /*
       * At this point the divisor is a true 'long' so we must use
       * Knuth's algorithm.
       *
       * Step D1: Normalize divisor and dividend (this makes our 'qhat'
       *        guesses more accurate):
       */
       for (shift=0; !SIGNBITON(divisor); shift++, divisor <<= 1) { ; }
       shift--;
       divisor >>= 1;
 
       if ((u1u2 >> (LONGSIZE - shift)) != 0 && shift != 0)
               abort("DLdiv:  dividend too large", 1);
       u1u2 = (u1u2 << shift) + ((shift == 0) ? 0 : u3u4 >> (LONGSIZE - shift));
       u3u4 <<= shift;
 
       /*
       * Step D2:  Begin Loop through digits, dividing u1,u2,u3 by v1,v2,
       *           then shifting U left by 1 digit:
       */
       v1 = HIGHDIGIT(divisor);
       v2 = LOWDIGIT(divisor);
       q3q4 = 0;
       u3 = HIGHDIGIT(u3u4);
 
       for (j=0; j < 2; j++) {
 
               /*
               * Step D3:  make a guess (qhat) at the next quotient denominator:
               */
               qhat = (HIGHDIGIT(u1u2) == v1) ? MAXSHORT : u1u2 / v1;
               /*
               * At this point Knuth would have us further refine our
               * guess, since we know qhat is too big if
               *
               *      v2 * qhat > ASSEMBLE(u1u2 % v, u3)
               *
               * That would make sense if u1u2 % v was easy to find, as it
               * would be in assembly language.  I ignore this step, and
               * repeat step D6 if qhat is too big.
               */
 
               /*
               * 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", 2);
               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); */
       return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void DLmult ( doublelong product,
ULONG  u,
ULONG  v 
)

Definition at line 110 of file arith.c.

{
  register ULONG u1, u2; /* the digits of u */
  register ULONG v1, v2; /* the digits of v */
  register unsigned int w1, w2, w3, w4; /* the digits of w */
  register ULONG t; /* temporary variable */
/* printf("DLmult(? ?, %x, %x)\n", u, v); */
  u1 = HIGHDIGIT(u);
  u2 = LOWDIGIT(u);
  v1 = HIGHDIGIT(v);
  v2 = LOWDIGIT(v);
 
  if (v2 == 0) w4 = w3 = w2 = 0;
  else
    {
    t = u2 * v2;
    w4 = LOWDIGIT(t);
    t = u1 * v2 + HIGHDIGIT(t);
    w3 = LOWDIGIT(t);
    w2 = HIGHDIGIT(t);
    }
 
  if (v1 == 0) w1 = 0;
  else
    {
    t = u2 * v1 + w3;
    w3 = LOWDIGIT(t);
    t = u1 * v1 + w2 + HIGHDIGIT(t);
    w2 = LOWDIGIT(t);
    w1 = HIGHDIGIT(t);
    }
 
  product->high = ASSEMBLE(w1, w2);
  product->low  = ASSEMBLE(w3, w4);
}

Here is the caller graph for this function:

void DLsub ( doublelong u,
doublelong v 
)

Definition at line 312 of file arith.c.

{
/* 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;
}
fractpel FPdiv ( fractpel  dividend,
fractpel  divisor 
)

Definition at line 372 of file arith.c.

{
       doublelong w;         /* result will be built here                    */
       int negative = FALSE; /* flag for sign bit                            */
       int maxshort = MAXSHORT; /* To avoid that overflow warning (RMz) */
 
       if (dividend < 0) {
               dividend = -dividend;
               negative = TRUE;
       }
       if (divisor < 0) {
               divisor = -divisor;
               negative = !negative;
       }
       w.low = dividend << FRACTBITS;
       w.high = dividend >> (LONGSIZE - FRACTBITS);
       DLdiv(&w, divisor);
       if (w.high != 0 || SIGNBITON(w.low)) {
               IfTrace2(TRUE,"FPdiv: overflow, %d/%d\n", dividend, divisor);
               w.low = TOFRACTPEL(maxshort);
       }
       return( (negative) ? -w.low : w.low);
}

Here is the call graph for this function:

fractpel FPmult ( fractpel  u,
fractpel  v 
)

Definition at line 340 of file arith.c.

{
  doublelong w;
  register int negative = FALSE; /* sign flag */
  int maxshort = MAXSHORT; /* To avoid that overflow warning (RMz) */
 
  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);
  if (w.high != 0 || SIGNBITON(w.low)) {
        IfTrace2(TRUE,"FPmult: overflow, %dx%d\n", u, v);
        w.low = TOFRACTPEL(maxshort);
  }
 
  return ((negative) ? -w.low : w.low);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Definition at line 406 of file arith.c.

{
       doublelong w;         /* result will be built here                    */
       int negative = FALSE;
       int maxshort = MAXSHORT; /* To avoid that overflow warning (RMz) */
       
       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);
       if (w.high != 0 || SIGNBITON(w.low)) {
               IfTrace3(TRUE,"FPstarslash: overflow, %d*%d/%d\n", a, b, c);
               w.low = TOFRACTPEL(maxshort);
       }
       return((negative) ? -w.low : w.low);
}

Here is the call graph for this function: