arith.c

远程桌面连接工具

/* $XConsortium: arith.c,v 1.4 94/03/22 19:08:54 gildea Exp $ */
/* Copyright International Business Machines, Corp. 1991
 * All Rights Reserved
 * Copyright Lexmark International, Inc. 1991
 * All Rights Reserved
 *
 * License to use, copy, modify, and distribute this software and its
 * documentation for any purpose and without fee is hereby granted,
 * provided that the above copyright notice appear in all copies and that
 * both that copyright notice and this permission notice appear in
 * supporting documentation, and that the name of IBM or Lexmark not be
 * used in advertising or publicity pertaining to distribution of the
 * software without specific, written prior permission.
 *
 * IBM AND LEXMARK PROVIDE THIS SOFTWARE "AS IS", WITHOUT ANY WARRANTIES OF
 * ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO ANY
 * IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE,
 * AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.  THE ENTIRE RISK AS TO THE
 * QUALITY AND PERFORMANCE OF THE SOFTWARE, INCLUDING ANY DUTY TO SUPPORT
 * OR MAINTAIN, BELONGS TO THE LICENSEE.  SHOULD ANY PORTION OF THE
 * SOFTWARE PROVE DEFECTIVE, THE LICENSEE (NOT IBM OR LEXMARK) ASSUMES THE
 * ENTIRE COST OF ALL SERVICING, REPAIR AND CORRECTION.  IN NO EVENT SHALL
 * IBM OR LEXMARK BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
 * DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
 * PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
 * THIS SOFTWARE.
 */
 /* ARITH    CWEB         V0006 ********                             */
/*
:h1.ARITH Module - Portable Module for Multiple Precision Fixed Point Arithmetic
 
This module provides division and multiplication of 64-bit fixed point
numbers.  (To be more precise, the module works on numbers that take
two 'longs' to store.  That is almost always equivalent to saying 64-bit
numbers.)
 
Note: it is frequently easy and desirable to recode these functions in
assembly language for the particular processor being used, because
assembly language, unlike C, will have 64-bit multiply products and
64-bit dividends.  This module is offered as a portable version.
 
&author. Jeffrey B. Lotspiech ([e-mail]) and Sten F. Andler
 
 
:h3.Include Files
 
The included files are:
*/
 
#include "objects.h"
#include "spaces.h"
#include "arith.h"
 
/*
:h3.
*/
/*SHARED LINE(S) ORIGINATED HERE*/
/*
Reference for all algorithms:  Donald E. Knuth, "The Art of Computer
Programming, Volume 2, Semi-Numerical Algorithms," Addison-Wesley Co.,
Massachusetts, 1969, pp. 229-279.
 
Knuth talks about a 'digit' being an arbitrary sized unit and a number
being a sequence of digits.  We'll take a digit to be a 'short'.
The following assumption must be valid for these algorithms to work:
:ol.
:li.A 'long' is two 'short's.
:eol.
The following code is INDEPENDENT of:
:ol.
:li.The actual size of a short.
:li.Whether shorts and longs are stored most significant byte
first or least significant byte first.
:eol.
 
SHORTSIZE is the number of bits in a short; LONGSIZE is the number of
bits in a long; MAXSHORT is the maximum unsigned short:
*/
/*SHARED LINE(S) ORIGINATED HERE*/
/*
ASSEMBLE concatenates two shorts to form a long:
*/
#define     ASSEMBLE(hi,lo)   ((((unsigned long)hi)<<SHORTSIZE)+(lo))
/*
HIGHDIGIT extracts the most significant short from a long; LOWDIGIT
extracts the least significant short from a long:
*/
#define     HIGHDIGIT(u)      ((u)>>SHORTSIZE)
#define     LOWDIGIT(u)       ((u)&MAXSHORT)
 
/*
SIGNBITON tests the high order bit of a long 'w':
*/
#define    SIGNBITON(w)   (((long)w)<0)
 
/*SHARED LINE(S) ORIGINATED HERE*/
 
/*
:h2.Double Long Arithmetic
 
:h3.DLmult() - Multiply Two Longs to Yield a Double Long
 
The two multiplicands must be positive.
*/
 
void DLmult(product, u, v)
  register doublelong *product;
  register unsigned long u;
  register unsigned long v;
{
#ifdef LONG64
/* printf("DLmult(? ?, %lx, %lx)\n", u, v); */
    *product = u*v;
/* printf("DLmult returns %lx\n", *product); */
#else
  register unsigned long u1, u2; /* the digits of u */
  register unsigned long v1, v2; /* the digits of v */
  register unsigned int w1, w2, w3, w4; /* the digits of w */
  register unsigned long 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);
#endif /* LONG64 else */
}
 
/*
:h2.DLdiv() - Divide Two Longs by One Long, Yielding Two Longs
 
Both the dividend and the divisor must be positive.
*/
 
void DLdiv(quotient, divisor)
       doublelong *quotient;       /* also where dividend is, originally     */
       unsigned long divisor;
{
#ifdef LONG64
/* printf("DLdiv(%lx %lx)\n", quotient, divisor); */
	*quotient /= divisor;
/* printf("DLdiv returns %lx\n", *quotient); */
#else
       register unsigned long u1u2 = quotient->high;
       register unsigned long 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 unsigned long 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");
       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.