usr/src/common/ficl/double.c

   1 /*
   2  * m a t h 6 4 . c
   3  * Forth Inspired Command Language - 64 bit math support routines
   4  * Authors: Michael A. Gauland (gaulandm@mdhost.cse.tek.com)
   5  *          Larry Hastings (larry@hastings.org)
   6  *          John Sadler (john_sadler@alum.mit.edu)
   7  * Created: 25 January 1998
   8  * Rev 2.03: Support for 128 bit DP math. This file really ouught to
   9  * be renamed!
  10  * $Id: double.c,v 1.2 2010/09/12 15:18:07 asau Exp $
  11  */
  12 /*
  13  * Copyright (c) 1997-2001 John Sadler (john_sadler@alum.mit.edu)
  14  * All rights reserved.
  15  *
  16  * Get the latest Ficl release at http://ficl.sourceforge.net
  17  *
  18  * I am interested in hearing from anyone who uses Ficl. If you have
  19  * a problem, a success story, a defect, an enhancement request, or
  20  * if you would like to contribute to the Ficl release, please
  21  * contact me by email at the address above.
  22  *
  23  * L I C E N S E  and  D I S C L A I M E R
  24  *
  25  * Redistribution and use in source and binary forms, with or without
  26  * modification, are permitted provided that the following conditions
  27  * are met:
  28  * 1. Redistributions of source code must retain the above copyright
  29  *    notice, this list of conditions and the following disclaimer.
  30  * 2. Redistributions in binary form must reproduce the above copyright
  31  *    notice, this list of conditions and the following disclaimer in the
  32  *    documentation and/or other materials provided with the distribution.
  33  *
  34  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  35  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  36  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  37  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  38  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  39  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  40  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  41  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  42  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  43  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  44  * SUCH DAMAGE.
  45  */
  46
  47 #include "ficl.h"
  48
  49 #if FICL_PLATFORM_HAS_2INTEGER
  50 ficl2UnsignedQR
  51 ficl2UnsignedDivide(ficl2Unsigned q, ficlUnsigned y)
  52 {
  53         ficl2UnsignedQR result;
  54
  55         result.quotient = q / y;
  56         /*
  57          * Once we have the quotient, it's cheaper to calculate the
  58          * remainder this way than with % (mod).  --lch
  59          */
  60         result.remainder  = (ficlInteger)(q - (result.quotient * y));
  61
  62         return (result);
  63 }
  64
  65 #else  /* FICL_PLATFORM_HAS_2INTEGER */
  66
  67 #define FICL_CELL_HIGH_BIT      ((uintmax_t)1 << (FICL_BITS_PER_CELL-1))
  68 #define UMOD_SHIFT              (FICL_BITS_PER_CELL / 2)
  69 #define UMOD_MASK               ((1L << (FICL_BITS_PER_CELL / 2)) - 1)
  70
  71 /*
  72  * ficl2IntegerIsNegative
  73  * Returns TRUE if the specified ficl2Unsigned has its sign bit set.
  74  */
  75 int
  76 ficl2IntegerIsNegative(ficl2Integer x)
  77 {
  78         return (x.high < 0);
  79 }
  80
  81 /*
  82  * ficl2IntegerNegate
  83  * Negates an ficl2Unsigned by complementing and incrementing.
  84  */
  85 ficl2Integer
  86 ficl2IntegerNegate(ficl2Integer x)
  87 {
  88         x.high = ~x.high;
  89         x.low = ~x.low;
  90         x.low ++;
  91         if (x.low == 0)
  92                 x.high++;
  93
  94         return (x);
  95 }
  96
  97 /*
  98  * ficl2UnsignedMultiplyAccumulate
  99  * Mixed precision multiply and accumulate primitive for number building.
 100  * Multiplies ficl2Unsigned u by ficlUnsigned mul and adds ficlUnsigned add.
 101  * Mul is typically the numeric base, and add represents a digit to be
 102  * appended to the growing number.
 103  * Returns the result of the operation
 104  */
 105 ficl2Unsigned
 106 ficl2UnsignedMultiplyAccumulate(ficl2Unsigned u, ficlUnsigned mul,
 107     ficlUnsigned add)
 108 {
 109         ficl2Unsigned resultLo = ficl2UnsignedMultiply(u.low, mul);
 110         ficl2Unsigned resultHi = ficl2UnsignedMultiply(u.high, mul);
 111         resultLo.high += resultHi.low;
 112         resultHi.low = resultLo.low + add;
 113
 114         if (resultHi.low < resultLo.low)
 115                 resultLo.high++;
 116
 117         resultLo.low = resultHi.low;
 118
 119         return (resultLo);
 120 }
 121
 122 /*
 123  * ficl2IntegerMultiply
 124  * Multiplies a pair of ficlIntegers and returns an ficl2Integer result.
 125  */
 126 ficl2Integer
 127 ficl2IntegerMultiply(ficlInteger x, ficlInteger y)
 128 {
 129         ficl2Unsigned prod;
 130         ficl2Integer result;
 131         int sign = 1;
 132
 133         if (x < 0) {
 134                 sign = -sign;
 135                 x = -x;
 136         }
 137
 138         if (y < 0) {
 139                 sign = -sign;
 140                 y = -y;
 141         }
 142
 143         prod = ficl2UnsignedMultiply(x, y);
 144         FICL_2INTEGER_SET(FICL_2UNSIGNED_GET_HIGH(prod),
 145             FICL_2UNSIGNED_GET_LOW(prod), result);
 146         if (sign > 0)
 147                 return (result);
 148         else
 149                 return (ficl2IntegerNegate(result));
 150 }
 151
 152 ficl2Integer
 153 ficl2IntegerDecrement(ficl2Integer x)
 154 {
 155         if (x.low == INTMAX_MIN)
 156                 x.high--;
 157         x.low--;
 158
 159         return (x);
 160 }
 161
 162 ficl2Unsigned
 163 ficl2UnsignedAdd(ficl2Unsigned x, ficl2Unsigned y)
 164 {
 165         ficl2Unsigned result;
 166
 167         result.high = x.high + y.high;
 168         result.low = x.low + y.low;
 169
 170         if (result.low < y.low)
 171                 result.high++;
 172
 173         return (result);
 174 }
 175
 176 /*
 177  * ficl2UnsignedMultiply
 178  * Contributed by:
 179  * Michael A. Gauland   gaulandm@mdhost.cse.tek.com
 180  */
 181 ficl2Unsigned
 182 ficl2UnsignedMultiply(ficlUnsigned x, ficlUnsigned y)
 183 {
 184         ficl2Unsigned result = { 0, 0 };
 185         ficl2Unsigned addend;
 186
 187         addend.low = y;
 188         addend.high = 0; /* No sign extension--arguments are unsigned */
 189
 190         while (x != 0) {
 191                 if (x & 1) {
 192                         result = ficl2UnsignedAdd(result, addend);
 193                 }
 194                 x >>= 1;
 195                 addend = ficl2UnsignedArithmeticShiftLeft(addend);
 196         }
 197         return (result);
 198 }
 199
 200 /*
 201  *                      ficl2UnsignedSubtract
 202  */
 203 ficl2Unsigned
 204 ficl2UnsignedSubtract(ficl2Unsigned x, ficl2Unsigned y)
 205 {
 206         ficl2Unsigned result;
 207
 208         result.high = x.high - y.high;
 209         result.low = x.low - y.low;
 210
 211         if (x.low < y.low) {
 212                 result.high--;
 213         }
 214
 215         return (result);
 216 }
 217
 218 /*
 219  * ficl2UnsignedArithmeticShiftLeft
 220  * 64 bit left shift
 221  */
 222 ficl2Unsigned
 223 ficl2UnsignedArithmeticShiftLeft(ficl2Unsigned x)
 224 {
 225         ficl2Unsigned result;
 226
 227         result.high = x.high << 1;
 228         if (x.low & FICL_CELL_HIGH_BIT) {
 229                 result.high++;
 230         }
 231
 232         result.low = x.low << 1;
 233
 234         return (result);
 235 }
 236
 237 /*
 238  * ficl2UnsignedArithmeticShiftRight
 239  * 64 bit right shift (unsigned - no sign extend)
 240  */
 241 ficl2Unsigned
 242 ficl2UnsignedArithmeticShiftRight(ficl2Unsigned x)
 243 {
 244         ficl2Unsigned result;
 245
 246         result.low = x.low >> 1;
 247         if (x.high & 1) {
 248                 result.low |= FICL_CELL_HIGH_BIT;
 249         }
 250
 251         result.high = x.high >> 1;
 252         return (result);
 253 }
 254
 255 /*
 256  * ficl2UnsignedOr
 257  * 64 bit bitwise OR
 258  */
 259 ficl2Unsigned
 260 ficl2UnsignedOr(ficl2Unsigned x, ficl2Unsigned y)
 261 {
 262         ficl2Unsigned result;
 263
 264         result.high = x.high | y.high;
 265         result.low = x.low | y.low;
 266
 267         return (result);
 268 }
 269
 270 /*
 271  * ficl2UnsignedCompare
 272  * Return -1 if x < y; 0 if x==y, and 1 if x > y.
 273  */
 274 int
 275 ficl2UnsignedCompare(ficl2Unsigned x, ficl2Unsigned y)
 276 {
 277         if (x.high > y.high)
 278                 return (1);
 279         if (x.high < y.high)
 280                 return (-1);
 281
 282         /* High parts are equal */
 283
 284         if (x.low > y.low)
 285                 return (1);
 286         else if (x.low < y.low)
 287                 return (-1);
 288
 289         return (0);
 290 }
 291
 292 /*
 293  * ficl2UnsignedDivide
 294  * Portable versions of ficl2Multiply and ficl2Divide in C
 295  * Contributed by:
 296  * Michael A. Gauland   gaulandm@mdhost.cse.tek.com
 297  */
 298 ficl2UnsignedQR
 299 ficl2UnsignedDivide(ficl2Unsigned q, ficlUnsigned y)
 300 {
 301         ficl2UnsignedQR result;
 302         ficl2Unsigned quotient;
 303         ficl2Unsigned subtrahend;
 304         ficl2Unsigned mask;
 305
 306         quotient.low = 0;
 307         quotient.high = 0;
 308
 309         subtrahend.low = y;
 310         subtrahend.high = 0;
 311
 312         mask.low = 1;
 313         mask.high = 0;
 314
 315         while ((ficl2UnsignedCompare(subtrahend, q) < 0) &&
 316             (subtrahend.high & FICL_CELL_HIGH_BIT) == 0) {
 317                 mask = ficl2UnsignedArithmeticShiftLeft(mask);
 318                 subtrahend = ficl2UnsignedArithmeticShiftLeft(subtrahend);
 319         }
 320
 321         while (mask.low != 0 || mask.high != 0) {
 322                 if (ficl2UnsignedCompare(subtrahend, q) <= 0) {
 323                         q = ficl2UnsignedSubtract(q, subtrahend);
 324                         quotient = ficl2UnsignedOr(quotient, mask);
 325                 }
 326                 mask = ficl2UnsignedArithmeticShiftRight(mask);
 327                 subtrahend = ficl2UnsignedArithmeticShiftRight(subtrahend);
 328         }
 329
 330         result.quotient = quotient;
 331         result.remainder = q.low;
 332         return (result);
 333 }
 334 #endif /* !FICL_PLATFORM_HAS_2INTEGER */
 335
 336 /*
 337  * ficl2IntegerDivideFloored
 338  *
 339  * FROM THE FORTH ANS...
 340  * Floored division is integer division in which the remainder carries
 341  * the sign of the divisor or is zero, and the quotient is rounded to
 342  * its arithmetic floor. Symmetric division is integer division in which
 343  * the remainder carries the sign of the dividend or is zero and the
 344  * quotient is the mathematical quotient rounded towards zero or
 345  * truncated. Examples of each are shown in tables 3.3 and 3.4.
 346  *
 347  * Table 3.3 - Floored Division Example
 348  * Dividend        Divisor Remainder       Quotient
 349  * --------        ------- ---------       --------
 350  *  10                7       3                1
 351  * -10                7       4               -2
 352  *  10               -7      -4               -2
 353  * -10               -7      -3                1
 354  *
 355  *
 356  * Table 3.4 - Symmetric Division Example
 357  * Dividend        Divisor Remainder       Quotient
 358  * --------        ------- ---------       --------
 359  *  10                7       3                1
 360  * -10                7      -3               -1
 361  *  10               -7       3               -1
 362  * -10               -7      -3                1
 363  */
 364 ficl2IntegerQR
 365 ficl2IntegerDivideFloored(ficl2Integer num, ficlInteger den)
 366 {
 367         ficl2IntegerQR qr;
 368         ficl2UnsignedQR uqr;
 369         ficl2Unsigned u;
 370         int signRem = 1;
 371         int signQuot = 1;
 372
 373         if (ficl2IntegerIsNegative(num)) {
 374                 num = ficl2IntegerNegate(num);
 375                 signQuot = -signQuot;
 376         }
 377
 378         if (den < 0) {
 379                 den = -den;
 380                 signRem = -signRem;
 381                 signQuot = -signQuot;
 382         }
 383
 384         FICL_2UNSIGNED_SET(FICL_2UNSIGNED_GET_HIGH(num),
 385             FICL_2UNSIGNED_GET_LOW(num), u);
 386         uqr = ficl2UnsignedDivide(u, (ficlUnsigned)den);
 387         qr = FICL_2UNSIGNEDQR_TO_2INTEGERQR(uqr);
 388         if (signQuot < 0) {
 389                 qr.quotient = ficl2IntegerNegate(qr.quotient);
 390                 if (qr.remainder != 0) {
 391                         qr.quotient = ficl2IntegerDecrement(qr.quotient);
 392                         qr.remainder = den - qr.remainder;
 393                 }
 394         }
 395
 396         if (signRem < 0)
 397                 qr.remainder = -qr.remainder;
 398
 399         return (qr);
 400 }
 401
 402 /*
 403  * ficl2IntegerDivideSymmetric
 404  * Divide an ficl2Unsigned by a ficlInteger and return a ficlInteger quotient
 405  * and a ficlInteger remainder. The absolute values of quotient and remainder
 406  * are not affected by the signs of the numerator and denominator
 407  * (the operation is symmetric on the number line)
 408  */
 409 ficl2IntegerQR
 410 ficl2IntegerDivideSymmetric(ficl2Integer num, ficlInteger den)
 411 {
 412         ficl2IntegerQR qr;
 413         ficl2UnsignedQR uqr;
 414         ficl2Unsigned u;
 415         int signRem = 1;
 416         int signQuot = 1;
 417
 418         if (ficl2IntegerIsNegative(num)) {
 419                 num = ficl2IntegerNegate(num);
 420                 signRem  = -signRem;
 421                 signQuot = -signQuot;
 422         }
 423
 424         if (den < 0) {
 425                 den = -den;
 426                 signQuot = -signQuot;
 427         }
 428
 429         FICL_2UNSIGNED_SET(FICL_2UNSIGNED_GET_HIGH(num),
 430             FICL_2UNSIGNED_GET_LOW(num), u);
 431         uqr = ficl2UnsignedDivide(u, (ficlUnsigned)den);
 432         qr = FICL_2UNSIGNEDQR_TO_2INTEGERQR(uqr);
 433         if (signRem < 0)
 434                 qr.remainder = -qr.remainder;
 435
 436         if (signQuot < 0)
 437                 qr.quotient = ficl2IntegerNegate(qr.quotient);
 438
 439         return (qr);
 440 }