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 <stdio.h>
00052 
00053 #include "types.h"
00054 #include "objects.h"
00055 #include "spaces.h"
00056 #include "arith.h"
00057  
00058 /*
00059 :h3.
00060 */
00061 /*SHARED LINE(S) ORIGINATED HERE*/
00062 /*
00063 Reference for all algorithms:  Donald E. Knuth, "The Art of Computer
00064 Programming, Volume 2, Semi-Numerical Algorithms," Addison-Wesley Co.,
00065 Massachusetts, 1969, pp. 229-279.
00066  
00067 Knuth talks about a 'digit' being an arbitrary sized unit and a number
00068 being a sequence of digits.  We'll take a digit to be a 'short'.
00069 The following assumption must be valid for these algorithms to work:
00070 :ol.
00071 :li.A 'long' is two 'short's.
00072 :eol.
00073 The following code is INDEPENDENT of:
00074 :ol.
00075 :li.The actual size of a short.
00076 :li.Whether shorts and longs are stored most significant byte
00077 first or least significant byte first.
00078 :eol.
00079  
00080 SHORTSIZE is the number of bits in a short; LONGSIZE is the number of
00081 bits in a long; MAXSHORT is the maximum unsigned short:
00082 */
00083 /*SHARED LINE(S) ORIGINATED HERE*/
00084 /*
00085 ASSEMBLE concatenates two shorts to form a long:
00086 */
00087 #define     ASSEMBLE(hi,lo)   ((((ULONG)hi)<<SHORTSIZE)+(lo))
00088 /*
00089 HIGHDIGIT extracts the most significant short from a long; LOWDIGIT
00090 extracts the least significant short from a long:
00091 */
00092 #define     HIGHDIGIT(u)      ((u)>>SHORTSIZE)
00093 #define     LOWDIGIT(u)       ((u)&MAXSHORT)
00094  
00095 /*
00096 SIGNBITON tests the high order bit of a long 'w':
00097 */
00098 #define    SIGNBITON(w)   (((LONG)w)<0)
00099  
00100 /*SHARED LINE(S) ORIGINATED HERE*/
00101  
00102 /*
00103 :h2.Double Long Arithmetic
00104  
00105 :h3.DLmult() - Multiply Two Longs to Yield a Double Long
00106  
00107 The two multiplicands must be positive.
00108 */
00109  
00110 void DLmult(product, u, v)
00111   register doublelong *product;
00112   register ULONG u;
00113   register ULONG v;
00114 {
00115   register ULONG u1, u2; /* the digits of u */
00116   register ULONG v1, v2; /* the digits of v */
00117   register unsigned int w1, w2, w3, w4; /* the digits of w */
00118   register ULONG t; /* temporary variable */
00119 /* printf("DLmult(? ?, %x, %x)\n", u, v); */
00120   u1 = HIGHDIGIT(u);
00121   u2 = LOWDIGIT(u);
00122   v1 = HIGHDIGIT(v);
00123   v2 = LOWDIGIT(v);
00124  
00125   if (v2 == 0) w4 = w3 = w2 = 0;
00126   else
00127     {
00128     t = u2 * v2;
00129     w4 = LOWDIGIT(t);
00130     t = u1 * v2 + HIGHDIGIT(t);
00131     w3 = LOWDIGIT(t);
00132     w2 = HIGHDIGIT(t);
00133     }
00134  
00135   if (v1 == 0) w1 = 0;
00136   else
00137     {
00138     t = u2 * v1 + w3;
00139     w3 = LOWDIGIT(t);
00140     t = u1 * v1 + w2 + HIGHDIGIT(t);
00141     w2 = LOWDIGIT(t);
00142     w1 = HIGHDIGIT(t);
00143     }
00144  
00145   product->high = ASSEMBLE(w1, w2);
00146   product->low  = ASSEMBLE(w3, w4);
00147 }
00148  
00149 /*
00150 :h2.DLdiv() - Divide Two Longs by One Long, Yielding Two Longs
00151  
00152 Both the dividend and the divisor must be positive.
00153 */
00154  
00155 void DLdiv(quotient, divisor)
00156        doublelong *quotient;       /* also where dividend is, originally     */
00157        ULONG divisor;
00158 {
00159        register ULONG u1u2 = quotient->high;
00160        register ULONG u3u4 = quotient->low;
00161        register LONG u3;     /* single digit of dividend                     */
00162        register int v1,v2;   /* divisor in registers                         */
00163        register LONG t;      /* signed copy of u1u2                          */
00164        register int qhat;    /* guess at the quotient digit                  */
00165        register ULONG q3q4;  /* low two digits of quotient           */
00166        register int shift;   /* holds the shift value for normalizing        */
00167        register int j;       /* loop variable                                */
00168  
00169 /* printf("DLdiv(%x %x, %x)\n", quotient->high, quotient->low, divisor); */
00170        /*
00171        * Knuth's algorithm works if the dividend is smaller than the
00172        * divisor.  We can get to that state quickly:
00173        */
00174        if (u1u2 >= divisor) {
00175                quotient->high = u1u2 / divisor;
00176                u1u2 %= divisor;
00177        }
00178        else
00179                quotient->high = 0;
00180  
00181        if (divisor <= MAXSHORT) {
00182  
00183                /*
00184                * This is the case where the divisor is contained in one
00185                * 'short'.  It is worthwhile making this fast:
00186                */
00187                u1u2 = ASSEMBLE(u1u2, HIGHDIGIT(u3u4));
00188                q3q4 = u1u2 / divisor;
00189                u1u2 %= divisor;
00190                u1u2 = ASSEMBLE(u1u2, LOWDIGIT(u3u4));
00191                quotient->low = ASSEMBLE(q3q4, u1u2 / divisor);
00192                return;
00193        }
00194  
00195  
00196        /*
00197        * At this point the divisor is a true 'long' so we must use
00198        * Knuth's algorithm.
00199        *
00200        * Step D1: Normalize divisor and dividend (this makes our 'qhat'
00201        *        guesses more accurate):
00202        */
00203        for (shift=0; !SIGNBITON(divisor); shift++, divisor <<= 1) { ; }
00204        shift--;
00205        divisor >>= 1;
00206  
00207        if ((u1u2 >> (LONGSIZE - shift)) != 0 && shift != 0)
00208                abort("DLdiv:  dividend too large", 1);
00209        u1u2 = (u1u2 << shift) + ((shift == 0) ? 0 : u3u4 >> (LONGSIZE - shift));
00210        u3u4 <<= shift;
00211  
00212        /*
00213        * Step D2:  Begin Loop through digits, dividing u1,u2,u3 by v1,v2,
00214        *           then shifting U left by 1 digit:
00215        */
00216        v1 = HIGHDIGIT(divisor);
00217        v2 = LOWDIGIT(divisor);
00218        q3q4 = 0;
00219        u3 = HIGHDIGIT(u3u4);
00220  
00221        for (j=0; j < 2; j++) {
00222  
00223                /*
00224                * Step D3:  make a guess (qhat) at the next quotient denominator:
00225                */
00226                qhat = (HIGHDIGIT(u1u2) == v1) ? MAXSHORT : u1u2 / v1;
00227                /*
00228                * At this point Knuth would have us further refine our
00229                * guess, since we know qhat is too big if
00230                *
00231                *      v2 * qhat > ASSEMBLE(u1u2 % v, u3)
00232                *
00233                * That would make sense if u1u2 % v was easy to find, as it
00234                * would be in assembly language.  I ignore this step, and
00235                * repeat step D6 if qhat is too big.
00236                */
00237  
00238                /*
00239                * Step D4: Multiply v1,v2 times qhat and subtract it from
00240                * u1,u2,u3:
00241                */
00242                u3 -= qhat * v2;
00243                /*
00244                * The high digit of u3 now contains the "borrow" for the
00245                * rest of the substraction from u1,u2.
00246                * Sometimes we can lose the sign bit with the above.
00247                * If so, we have to force the high digit negative:
00248                */
00249                t = HIGHDIGIT(u3);
00250                if (t > 0)
00251                        t |= -1 << SHORTSIZE;
00252                t += u1u2 - qhat * v1;
00253 /* printf("..>divide step qhat=%x t=%x u3=%x u1u2=%x v1=%x v2=%x\n",
00254                              qhat, t, u3, u1u2, v1, v2); */
00255                while (t < 0) {  /* Test is Step D5.                          */
00256  
00257                        /*
00258                        * D6: Oops, qhat was too big.  Add back in v1,v2 and
00259                        * decrease qhat by 1:
00260                        */
00261                        u3 = LOWDIGIT(u3) + v2;
00262                        t += HIGHDIGIT(u3) + v1;
00263                        qhat--;
00264 /* printf("..>>qhat correction t=%x u3=%x qhat=%x\n", t, u3, qhat); */
00265                }
00266                /*
00267                * Step D7:  shift U left one digit and loop:
00268                */
00269                u1u2 = t;
00270                if (HIGHDIGIT(u1u2) != 0)
00271                        abort("divide algorithm error", 2);
00272                u1u2 = ASSEMBLE(u1u2, LOWDIGIT(u3));
00273                u3 = LOWDIGIT(u3u4);
00274                q3q4 = ASSEMBLE(q3q4, qhat);
00275        }
00276        quotient->low = q3q4;
00277 /* printf("DLdiv returns %x %x\n", quotient->high, quotient->low); */
00278        return;
00279 }
00280  
00281 /*
00282 :h3.DLadd() - Add Two Double Longs
00283  
00284 In this case, the doublelongs may be signed.  The algorithm takes the
00285 piecewise sum of the high and low longs, with the possibility that the
00286 high should be incremented if there is a carry out of the low.  How to
00287 tell if there is a carry?  Alex Harbury suggested that if the sum of
00288 the lows is less than the max of the lows, there must have been a
00289 carry.  Conversely, if there was a carry, the sum of the lows must be
00290 less than the max of the lows.  So, the test is "if and only if".
00291 */
00292  
00293 void DLadd(u, v)
00294        doublelong *u;        /* u = u + v                                    */
00295        doublelong *v;
00296 {
00297        register ULONG lowmax = TYPE1_MAX(u->low, v->low);
00298  
00299 /* printf("DLadd(%x %x, %x %x)\n", u->high, u->low, v->high, v->low); */
00300        u->high += v->high;
00301        u->low += v->low;
00302        if (lowmax > u->low)
00303                u->high++;
00304 }
00305 /*
00306 :h3.DLsub() - Subtract Two Double Longs
00307  
00308 Testing for a borrow is even easier.  If the v.low is greater than
00309 u.low, there must be a borrow.
00310 */
00311  
00312 void DLsub(u, v)
00313        doublelong *u;        /* u = u - v                                    */
00314        doublelong *v;
00315 {
00316 /* printf("DLsub(%x %x, %x %x)\n", u->high, u->low, v->high, v->low);*/
00317        u->high -= v->high;
00318        if (v->low > u->low)
00319                u->high--;
00320        u->low -= v->low;
00321 }
00322 /*
00323 :h3.DLrightshift() - Macro to Shift Double Long Right by N
00324 */
00325  
00326 /*SHARED LINE(S) ORIGINATED HERE*/
00327  
00328 /*
00329 :h2.Fractional Pel Arithmetic
00330 */
00331 /*
00332 :h3.FPmult() - Multiply Two Fractional Pel Values
00333  
00334 This funtion first calculates w = u * v to "doublelong" precision.
00335 It then shifts w right by FRACTBITS bits, and checks that no
00336 overflow will occur when the resulting value is passed back as
00337 a fractpel.
00338 */
00339  
00340 fractpel FPmult(u, v)
00341   register fractpel u,v;
00342 {
00343   doublelong w;
00344   register int negative = FALSE; /* sign flag */
00345   int maxshort = MAXSHORT; /* To avoid that overflow warning (RMz) */
00346  
00347   if ((u == 0) || (v == 0)) return (0);
00348  
00349  
00350   if (u < 0) {u = -u; negative = TRUE;}
00351   if (v < 0) {v = -v; negative = !negative;}
00352  
00353   if (u == TOFRACTPEL(1)) return ((negative) ? -v : v);
00354   if (v == TOFRACTPEL(1)) return ((negative) ? -u : u);
00355  
00356   DLmult(&w, u, v);
00357   DLrightshift(w, FRACTBITS);
00358   if (w.high != 0 || SIGNBITON(w.low)) {
00359         IfTrace2(TRUE,"FPmult: overflow, %dx%d\n", u, v);
00360         w.low = TOFRACTPEL(maxshort);
00361   }
00362  
00363   return ((negative) ? -w.low : w.low);
00364 }
00365  
00366 /*
00367 :h3.FPdiv() - Divide Two Fractional Pel Values
00368  
00369 These values may be signed.  The function returns the quotient.
00370 */
00371  
00372 fractpel FPdiv(dividend, divisor)
00373        register fractpel dividend;
00374        register fractpel divisor;
00375 {
00376        doublelong w;         /* result will be built here                    */
00377        int negative = FALSE; /* flag for sign bit                            */
00378        int maxshort = MAXSHORT; /* To avoid that overflow warning (RMz) */
00379  
00380        if (dividend < 0) {
00381                dividend = -dividend;
00382                negative = TRUE;
00383        }
00384        if (divisor < 0) {
00385                divisor = -divisor;
00386                negative = !negative;
00387        }
00388        w.low = dividend << FRACTBITS;
00389        w.high = dividend >> (LONGSIZE - FRACTBITS);
00390        DLdiv(&w, divisor);
00391        if (w.high != 0 || SIGNBITON(w.low)) {
00392                IfTrace2(TRUE,"FPdiv: overflow, %d/%d\n", dividend, divisor);
00393                w.low = TOFRACTPEL(maxshort);
00394        }
00395        return( (negative) ? -w.low : w.low);
00396 }
00397  
00398 /*
00399 :h3.FPstarslash() - Multiply then Divide
00400  
00401 Borrowing a chapter from the language Forth, it is useful to define
00402 an operator that first multiplies by one constant then divides by
00403 another, keeping the intermediate result in extended precision.
00404 */
00405  
00406 fractpel FPstarslash(a, b, c)
00407        register fractpel a,b,c;  /* result = a * b / c                       */
00408 {
00409        doublelong w;         /* result will be built here                    */
00410        int negative = FALSE;
00411        int maxshort = MAXSHORT; /* To avoid that overflow warning (RMz) */
00412        
00413        if (a < 0) { a = -a; negative = TRUE; }
00414        if (b < 0) { b = -b; negative = !negative; }
00415        if (c < 0) { c = -c; negative = !negative; }
00416  
00417        DLmult(&w, a, b);
00418        DLdiv(&w, c);
00419        if (w.high != 0 || SIGNBITON(w.low)) {
00420                IfTrace3(TRUE,"FPstarslash: overflow, %d*%d/%d\n", a, b, c);
00421                w.low = TOFRACTPEL(maxshort);
00422        }
00423        return((negative) ? -w.low : w.low);
00424 }