Back to index

tetex-bin  3.0
arith.c
Go to the documentation of this file.
00001 /* $XConsortium: arith.c,v 1.2 91/10/10 11:14:06 rws Exp $ */
00002 /* Copyright International Business Machines, Corp. 1991
00003  * All Rights Reserved
00004  * Copyright Lexmark International, Inc. 1991
00005  * All Rights Reserved
00006  *
00007  * License to use, copy, modify, and distribute this software and its
00008  * documentation for any purpose and without fee is hereby granted,
00009  * provided that the above copyright notice appear in all copies and that
00010  * both that copyright notice and this permission notice appear in
00011  * supporting documentation, and that the name of IBM or Lexmark not be
00012  * used in advertising or publicity pertaining to distribution of the
00013  * software without specific, written prior permission.
00014  *
00015  * IBM AND LEXMARK PROVIDE THIS SOFTWARE "AS IS", WITHOUT ANY WARRANTIES OF
00016  * ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO ANY
00017  * IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE,
00018  * AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.  THE ENTIRE RISK AS TO THE
00019  * QUALITY AND PERFORMANCE OF THE SOFTWARE, INCLUDING ANY DUTY TO SUPPORT
00020  * OR MAINTAIN, BELONGS TO THE LICENSEE.  SHOULD ANY PORTION OF THE
00021  * SOFTWARE PROVE DEFECTIVE, THE LICENSEE (NOT IBM OR LEXMARK) ASSUMES THE
00022  * ENTIRE COST OF ALL SERVICING, REPAIR AND CORRECTION.  IN NO EVENT SHALL
00023  * IBM OR LEXMARK BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
00024  * DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
00025  * PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
00026  * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
00027  * THIS SOFTWARE.
00028  */
00029  /* ARITH    CWEB         V0006 ********                             */
00030 /*
00031 :h1.ARITH Module - Portable Module for Multiple Precision Fixed Point Arithmetic
00032  
00033 This module provides division and multiplication of 64-bit fixed point
00034 numbers.  (To be more precise, the module works on numbers that take
00035 two 'longs' to store.  That is almost always equivalent to saying 64-bit
00036 numbers.)
00037  
00038 Note: it is frequently easy and desirable to recode these functions in
00039 assembly language for the particular processor being used, because
00040 assembly language, unlike C, will have 64-bit multiply products and
00041 64-bit dividends.  This module is offered as a portable version.
00042  
00043 &author. Jeffrey B. Lotspiech (lotspiech@almaden.ibm.com) and Sten F. Andler
00044  
00045  
00046 :h3.Include Files
00047  
00048 The included files are:
00049 */
00050  
00051 #include "types.h"
00052 #include "objects.h"
00053 #include "spaces.h"
00054 #include "arith.h"
00055  
00056 /*
00057 :h3.
00058 */
00059 /*SHARED LINE(S) ORIGINATED HERE*/
00060 /*
00061 Reference for all algorithms:  Donald E. Knuth, "The Art of Computer
00062 Programming, Volume 2, Semi-Numerical Algorithms," Addison-Wesley Co.,
00063 Massachusetts, 1969, pp. 229-279.
00064  
00065 Knuth talks about a 'digit' being an arbitrary sized unit and a number
00066 being a sequence of digits.  We'll take a digit to be a 'short'.
00067 The following assumption must be valid for these algorithms to work:
00068 :ol.
00069 :li.A 'long' is two 'short's.
00070 :eol.
00071 The following code is INDEPENDENT of:
00072 :ol.
00073 :li.The actual size of a short.
00074 :li.Whether shorts and longs are stored most significant byte
00075 first or least significant byte first.
00076 :eol.
00077  
00078 SHORTSIZE is the number of bits in a short; LONGSIZE is the number of
00079 bits in a long; MAXSHORT is the maximum unsigned short:
00080 */
00081 /*SHARED LINE(S) ORIGINATED HERE*/
00082 /*
00083 ASSEMBLE concatenates two shorts to form a long:
00084 */
00085 #define     ASSEMBLE(hi,lo)   ((((ULONG)hi)<<SHORTSIZE)+(lo))
00086 /*
00087 HIGHDIGIT extracts the most significant short from a long; LOWDIGIT
00088 extracts the least significant short from a long:
00089 */
00090 #define     HIGHDIGIT(u)      ((u)>>SHORTSIZE)
00091 #define     LOWDIGIT(u)       ((u)&MAXSHORT)
00092  
00093 /*
00094 SIGNBITON tests the high order bit of a long 'w':
00095 */
00096 #define    SIGNBITON(w)   (((LONG)w)<0)
00097  
00098 /*SHARED LINE(S) ORIGINATED HERE*/
00099  
00100 /*
00101 :h2.Double Long Arithmetic
00102  
00103 :h3.DLmult() - Multiply Two Longs to Yield a Double Long
00104  
00105 The two multiplicands must be positive.
00106 */
00107  
00108 void DLmult(product, u, v)
00109   register doublelong *product;
00110   register ULONG u;
00111   register ULONG v;
00112 {
00113   register ULONG u1, u2; /* the digits of u */
00114   register ULONG v1, v2; /* the digits of v */
00115   register unsigned int w1, w2, w3, w4; /* the digits of w */
00116   register ULONG t; /* temporary variable */
00117 /* printf("DLmult(? ?, %x, %x)\n", u, v); */
00118   u1 = HIGHDIGIT(u);
00119   u2 = LOWDIGIT(u);
00120   v1 = HIGHDIGIT(v);
00121   v2 = LOWDIGIT(v);
00122  
00123   if (v2 == 0) w4 = w3 = w2 = 0;
00124   else
00125     {
00126     t = u2 * v2;
00127     w4 = LOWDIGIT(t);
00128     t = u1 * v2 + HIGHDIGIT(t);
00129     w3 = LOWDIGIT(t);
00130     w2 = HIGHDIGIT(t);
00131     }
00132  
00133   if (v1 == 0) w1 = 0;
00134   else
00135     {
00136     t = u2 * v1 + w3;
00137     w3 = LOWDIGIT(t);
00138     t = u1 * v1 + w2 + HIGHDIGIT(t);
00139     w2 = LOWDIGIT(t);
00140     w1 = HIGHDIGIT(t);
00141     }
00142  
00143   product->high = ASSEMBLE(w1, w2);
00144   product->low  = ASSEMBLE(w3, w4);
00145 }
00146  
00147 /*
00148 :h2.DLdiv() - Divide Two Longs by One Long, Yielding Two Longs
00149  
00150 Both the dividend and the divisor must be positive.
00151 */
00152  
00153 void DLdiv(quotient, divisor)
00154        doublelong *quotient;       /* also where dividend is, originally     */
00155        ULONG divisor;
00156 {
00157        register ULONG u1u2 = quotient->high;
00158        register ULONG u3u4 = quotient->low;
00159        register LONG u3;     /* single digit of dividend                     */
00160        register int v1,v2;   /* divisor in registers                         */
00161        register LONG t;      /* signed copy of u1u2                          */
00162        register int qhat;    /* guess at the quotient digit                  */
00163        register ULONG q3q4;  /* low two digits of quotient           */
00164        register int shift;   /* holds the shift value for normalizing        */
00165        register int j;       /* loop variable                                */
00166  
00167 /* printf("DLdiv(%x %x, %x)\n", quotient->high, quotient->low, divisor); */
00168        /*
00169        * Knuth's algorithm works if the dividend is smaller than the
00170        * divisor.  We can get to that state quickly:
00171        */
00172        if (u1u2 >= divisor) {
00173                quotient->high = u1u2 / divisor;
00174                u1u2 %= divisor;
00175        }
00176        else
00177                quotient->high = 0;
00178  
00179        if (divisor <= MAXSHORT) {
00180  
00181                /*
00182                * This is the case where the divisor is contained in one
00183                * 'short'.  It is worthwhile making this fast:
00184                */
00185                u1u2 = ASSEMBLE(u1u2, HIGHDIGIT(u3u4));
00186                q3q4 = u1u2 / divisor;
00187                u1u2 %= divisor;
00188                u1u2 = ASSEMBLE(u1u2, LOWDIGIT(u3u4));
00189                quotient->low = ASSEMBLE(q3q4, u1u2 / divisor);
00190                return;
00191        }
00192  
00193  
00194        /*
00195        * At this point the divisor is a true 'long' so we must use
00196        * Knuth's algorithm.
00197        *
00198        * Step D1: Normalize divisor and dividend (this makes our 'qhat'
00199        *        guesses more accurate):
00200        */
00201        for (shift=0; !SIGNBITON(divisor); shift++, divisor <<= 1) { ; }
00202        shift--;
00203        divisor >>= 1;
00204  
00205        if ((u1u2 >> (LONGSIZE - shift)) != 0 && shift != 0)
00206                t1_abort("DLdiv:  dividend too large");
00207        u1u2 = (u1u2 << shift) + ((shift == 0) ? 0 : u3u4 >> (LONGSIZE - shift));
00208        u3u4 <<= shift;
00209  
00210        /*
00211        * Step D2:  Begin Loop through digits, dividing u1,u2,u3 by v1,v2,
00212        *           then shifting U left by 1 digit:
00213        */
00214        v1 = HIGHDIGIT(divisor);
00215        v2 = LOWDIGIT(divisor);
00216        q3q4 = 0;
00217        u3 = HIGHDIGIT(u3u4);
00218  
00219        for (j=0; j < 2; j++) {
00220  
00221                /*
00222                * Step D3:  make a guess (qhat) at the next quotient denominator:
00223                */
00224                qhat = (HIGHDIGIT(u1u2) == v1) ? MAXSHORT : u1u2 / v1;
00225                /*
00226                * At this point Knuth would have us further refine our
00227                * guess, since we know qhat is too big if
00228                *
00229                *      v2 * qhat > ASSEMBLE(u1u2 % v, u3)
00230                *
00231                * That would make sense if u1u2 % v was easy to find, as it
00232                * would be in assembly language.  I ignore this step, and
00233                * repeat step D6 if qhat is too big.
00234                */
00235  
00236                /*
00237                * Step D4: Multiply v1,v2 times qhat and subtract it from
00238                * u1,u2,u3:
00239                */
00240                u3 -= qhat * v2;
00241                /*
00242                * The high digit of u3 now contains the "borrow" for the
00243                * rest of the substraction from u1,u2.
00244                * Sometimes we can lose the sign bit with the above.
00245                * If so, we have to force the high digit negative:
00246                */
00247                t = HIGHDIGIT(u3);
00248                if (t > 0)
00249                        t |= -1 << SHORTSIZE;
00250                t += u1u2 - qhat * v1;
00251 /* printf("..>divide step qhat=%x t=%x u3=%x u1u2=%x v1=%x v2=%x\n",
00252                              qhat, t, u3, u1u2, v1, v2); */
00253                while (t < 0) {  /* Test is Step D5.                          */
00254  
00255                        /*
00256                        * D6: Oops, qhat was too big.  Add back in v1,v2 and
00257                        * decrease qhat by 1:
00258                        */
00259                        u3 = LOWDIGIT(u3) + v2;
00260                        t += HIGHDIGIT(u3) + v1;
00261                        qhat--;
00262 /* printf("..>>qhat correction t=%x u3=%x qhat=%x\n", t, u3, qhat); */
00263                }
00264                /*
00265                * Step D7:  shift U left one digit and loop:
00266                */
00267                u1u2 = t;
00268                if (HIGHDIGIT(u1u2) != 0)
00269                        t1_abort("divide algorithm error");
00270                u1u2 = ASSEMBLE(u1u2, LOWDIGIT(u3));
00271                u3 = LOWDIGIT(u3u4);
00272                q3q4 = ASSEMBLE(q3q4, qhat);
00273        }
00274        quotient->low = q3q4;
00275 /* printf("DLdiv returns %x %x\n", quotient->high, quotient->low); */
00276        return;
00277 }
00278  
00279 /*
00280 :h3.DLadd() - Add Two Double Longs
00281  
00282 In this case, the doublelongs may be signed.  The algorithm takes the
00283 piecewise sum of the high and low longs, with the possibility that the
00284 high should be incremented if there is a carry out of the low.  How to
00285 tell if there is a carry?  Alex Harbury suggested that if the sum of
00286 the lows is less than the max of the lows, there must have been a
00287 carry.  Conversely, if there was a carry, the sum of the lows must be
00288 less than the max of the lows.  So, the test is "if and only if".
00289 */
00290  
00291 void DLadd(u, v)
00292        doublelong *u;        /* u = u + v                                    */
00293        doublelong *v;
00294 {
00295        register ULONG lowmax = MAX(u->low, v->low);
00296  
00297 /* printf("DLadd(%x %x, %x %x)\n", u->high, u->low, v->high, v->low); */
00298        u->high += v->high;
00299        u->low += v->low;
00300        if (lowmax > u->low)
00301                u->high++;
00302 }
00303 /*
00304 :h3.DLsub() - Subtract Two Double Longs
00305  
00306 Testing for a borrow is even easier.  If the v.low is greater than
00307 u.low, there must be a borrow.
00308 */
00309  
00310 void DLsub(u, v)
00311        doublelong *u;        /* u = u - v                                    */
00312        doublelong *v;
00313 {
00314 /* printf("DLsub(%x %x, %x %x)\n", u->high, u->low, v->high, v->low);*/
00315        u->high -= v->high;
00316        if (v->low > u->low)
00317                u->high--;
00318        u->low -= v->low;
00319 }
00320 /*
00321 :h3.DLrightshift() - Macro to Shift Double Long Right by N
00322 */
00323  
00324 /*SHARED LINE(S) ORIGINATED HERE*/
00325  
00326 /*
00327 :h2.Fractional Pel Arithmetic
00328 */
00329 /*
00330 :h3.FPmult() - Multiply Two Fractional Pel Values
00331  
00332 This funtion first calculates w = u * v to "doublelong" precision.
00333 It then shifts w right by FRACTBITS bits, and checks that no
00334 overflow will occur when the resulting value is passed back as
00335 a fractpel.
00336 */
00337  
00338 fractpel FPmult(u, v)
00339   register fractpel u,v;
00340 {
00341   doublelong w;
00342   register int negative = FALSE; /* sign flag */
00343  
00344   if ((u == 0) || (v == 0)) return (0);
00345  
00346  
00347   if (u < 0) {u = -u; negative = TRUE;}
00348   if (v < 0) {v = -v; negative = !negative;}
00349  
00350   if (u == TOFRACTPEL(1)) return ((negative) ? -v : v);
00351   if (v == TOFRACTPEL(1)) return ((negative) ? -u : u);
00352  
00353   DLmult(&w, u, v);
00354   DLrightshift(w, FRACTBITS);
00355   if (w.high != 0 || SIGNBITON(w.low)) {
00356         IfTrace2(TRUE,"FPmult: overflow, %dlx%dl\n", u, v);
00357         w.low = TOFRACTPEL(MAXSHORT);
00358   }
00359  
00360   return ((negative) ? -w.low : w.low);
00361 }
00362  
00363 /*
00364 :h3.FPdiv() - Divide Two Fractional Pel Values
00365  
00366 These values may be signed.  The function returns the quotient.
00367 */
00368  
00369 fractpel FPdiv(dividend, divisor)
00370        register fractpel dividend;
00371        register fractpel divisor;
00372 {
00373        doublelong w;         /* result will be built here                    */
00374        int negative = FALSE; /* flag for sign bit                            */
00375  
00376        if (dividend < 0) {
00377                dividend = -dividend;
00378                negative = TRUE;
00379        }
00380        if (divisor < 0) {
00381                divisor = -divisor;
00382                negative = !negative;
00383        }
00384        w.low = dividend << FRACTBITS;
00385        w.high = dividend >> (LONGSIZE - FRACTBITS);
00386        DLdiv(&w, divisor);
00387        if (w.high != 0 || SIGNBITON(w.low)) {
00388                IfTrace2(TRUE,"FPdiv: overflow, %dl/%dl\n", dividend, divisor);
00389                w.low = TOFRACTPEL(MAXSHORT);
00390        }
00391        return( (negative) ? -w.low : w.low);
00392 }
00393  
00394 /*
00395 :h3.FPstarslash() - Multiply then Divide
00396  
00397 Borrowing a chapter from the language Forth, it is useful to define
00398 an operator that first multiplies by one constant then divides by
00399 another, keeping the intermediate result in extended precision.
00400 */
00401  
00402 fractpel FPstarslash(a, b, c)
00403        register fractpel a,b,c;  /* result = a * b / c                       */
00404 {
00405        doublelong w;         /* result will be built here                    */
00406        int negative = FALSE;
00407  
00408        if (a < 0) { a = -a; negative = TRUE; }
00409        if (b < 0) { b = -b; negative = !negative; }
00410        if (c < 0) { c = -c; negative = !negative; }
00411  
00412        DLmult(&w, a, b);
00413        DLdiv(&w, c);
00414        if (w.high != 0 || SIGNBITON(w.low)) {
00415                IfTrace3(TRUE,"FPstarslash: overflow, %dl*%dl/%dl\n", a, b, c);
00416                w.low = TOFRACTPEL(MAXSHORT);
00417        }
00418        return((negative) ? -w.low : w.low);
00419 }