/* $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 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 (lotspiech@almaden.ibm.com) and Sten F. Andler * * * * 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: * A 'long' is two 'short's. * * The following code is INDEPENDENT of: * The actual size of a short. * Whether shorts and longs are stored most significant byte * first or least significant byte first. */ /* * Include Files * * The included files are: */ #ifndef T1GST #include "global.h" #endif #include /** * Structure Definitions **/ typedef struct { long high; unsigned long low; } doublelong; /* * SHORTSIZE is the number of bits in a short; LONGSIZE is the number of * bits in a long; MAXSHORT is the maximum unsigned short: */ #define SHORTSIZE (sizeof(short)*8) #define LONGSIZE (SHORTSIZE*2) #define MAXSHORT ((1<>SHORTSIZE) #define LOWDIGIT(u) ((u)&MAXSHORT) /* * SIGNBITON tests the high order bit of a long 'w': */ #define SIGNBITON(w) (((long)w)<0) /** * Local prototypes **/ static void DLmult(doublelong *product, unsigned long u, unsigned long v); /* * DLrightshift() - Macro to Shift Double Long Right by N */ #define DLrightshift(dl,N) { \ dl.low = (dl.low >> N) + (((unsigned long) dl.high) << (LONGSIZE - N)); \ dl.high >>= N; \ } /* * Double Long Arithmetic * * DLmult() - Multiply Two Longs to Yield a Double Long * * The two multiplicands must be positive. */ // //OPTIMIZATION BY AMISH 4/26/93 //In my test cases, this function was always called with //u OR v == 1 (NOT TOFRACTPEL(1), but rather a 32 bit # //with the first bit on. - (important distinction) // //So, calling DLmult is silly in these cases, since we //can just set product.low = u or v (whichever is not 1) //and product.high = 0. // //NOTE: This resulted in a barely measurable speedup // static void DLmult(doublelong *product, unsigned long u, unsigned long v) { unsigned long u1, u2; /* the digits of u */ unsigned long v1, v2; /* the digits of v */ unsigned int w1, w2, w3, w4; /* the digits of w */ unsigned long t; /* temporary variable */ //THE FOLLOWING TWO IF STATEMENTS ADDED BY AMISH 4/26/93 if (u == 1) { product->high = 0; product->low = v; return; } if (v == 1) { product->high = 0; product->low = u; return; } 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); } /* * Fractional Pel Arithmetic */ /* * 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(fractpel u, fractpel v) { doublelong w; int negative = FALSE; /* sign flag */ 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, %px%p\n", u, v); w.low = TOFRACTPEL(MAXSHORT); } return (fractpel)((negative) ? -w.low : w.low); }