sys/boot/ficl/math64.c

   1 /*******************************************************************
   2 ** m a t h 6 4 . c
   3 ** Forth Inspired Command Language - 64 bit math support routines
   4 ** Author: John Sadler (john_sadler@alum.mit.edu)
   5 ** Created: 25 January 1998
   6 ** Rev 2.03: Support for 128 bit DP math. This file really ouught to
   7 ** be renamed!
   8 ** $Id: math64.c,v 1.9 2001/12/05 07:21:34 jsadler Exp $
   9 *******************************************************************/
  10 /*
  11 ** Copyright (c) 1997-2001 John Sadler (john_sadler@alum.mit.edu)
  12 ** All rights reserved.
  13 **
  14 ** Get the latest Ficl release at http://ficl.sourceforge.net
  15 **
  16 ** I am interested in hearing from anyone who uses ficl. If you have
  17 ** a problem, a success story, a defect, an enhancement request, or
  18 ** if you would like to contribute to the ficl release, please
  19 ** contact me by email at the address above.
  20 **
  21 ** L I C E N S E  and  D I S C L A I M E R
  22 **
  23 ** Redistribution and use in source and binary forms, with or without
  24 ** modification, are permitted provided that the following conditions
  25 ** are met:
  26 ** 1. Redistributions of source code must retain the above copyright
  27 **    notice, this list of conditions and the following disclaimer.
  28 ** 2. Redistributions in binary form must reproduce the above copyright
  29 **    notice, this list of conditions and the following disclaimer in the
  30 **    documentation and/or other materials provided with the distribution.
  31 **
  32 ** THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  33 ** ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  34 ** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  35 ** ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  36 ** FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  37 ** DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  38 ** OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  39 ** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  40 ** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  41 ** OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  42 ** SUCH DAMAGE.
  43 */
  44
  45 /*
  46  * $FreeBSD: src/sys/boot/ficl/math64.c,v 1.4 2002/04/09 17:45:11 dcs Exp $
  47  * $DragonFly: src/sys/boot/ficl/math64.c,v 1.3 2003/11/10 06:08:33 dillon Exp $
  48  */
  49
  50 #include "ficl.h"
  51 #include "math64.h"
  52
  53
  54 /**************************************************************************
  55                         m 6 4 A b s
  56 ** Returns the absolute value of an DPINT
  57 **************************************************************************/
  58 DPINT m64Abs(DPINT x)
  59 {
  60     if (m64IsNegative(x))
  61         x = m64Negate(x);
  62
  63     return x;
  64 }
  65
  66
  67 /**************************************************************************
  68                         m 6 4 F l o o r e d D i v I
  69 **
  70 ** FROM THE FORTH ANS...
  71 ** Floored division is integer division in which the remainder carries
  72 ** the sign of the divisor or is zero, and the quotient is rounded to
  73 ** its arithmetic floor. Symmetric division is integer division in which
  74 ** the remainder carries the sign of the dividend or is zero and the
  75 ** quotient is the mathematical quotient rounded towards zero or
  76 ** truncated. Examples of each are shown in tables 3.3 and 3.4.
  77 **
  78 ** Table 3.3 - Floored Division Example
  79 ** Dividend        Divisor Remainder       Quotient
  80 ** --------        ------- ---------       --------
  81 **  10                7       3                1
  82 ** -10                7       4               -2
  83 **  10               -7      -4               -2
  84 ** -10               -7      -3                1
  85 **
  86 **
  87 ** Table 3.4 - Symmetric Division Example
  88 ** Dividend        Divisor Remainder       Quotient
  89 ** --------        ------- ---------       --------
  90 **  10                7       3                1
  91 ** -10                7      -3               -1
  92 **  10               -7       3               -1
  93 ** -10               -7      -3                1
  94 **************************************************************************/
  95 INTQR m64FlooredDivI(DPINT num, FICL_INT den)
  96 {
  97     INTQR qr;
  98     UNSQR uqr;
  99     int signRem = 1;
 100     int signQuot = 1;
 101
 102     if (m64IsNegative(num))
 103     {
 104         num = m64Negate(num);
 105         signQuot = -signQuot;
 106     }
 107
 108     if (den < 0)
 109     {
 110         den      = -den;
 111         signRem  = -signRem;
 112         signQuot = -signQuot;
 113     }
 114
 115     uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
 116     qr = m64CastQRUI(uqr);
 117     if (signQuot < 0)
 118     {
 119         qr.quot = -qr.quot;
 120         if (qr.rem != 0)
 121         {
 122             qr.quot--;
 123             qr.rem = den - qr.rem;
 124         }
 125     }
 126
 127     if (signRem < 0)
 128         qr.rem = -qr.rem;
 129
 130     return qr;
 131 }
 132
 133
 134 /**************************************************************************
 135                         m 6 4 I s N e g a t i v e
 136 ** Returns TRUE if the specified DPINT has its sign bit set.
 137 **************************************************************************/
 138 int m64IsNegative(DPINT x)
 139 {
 140     return (x.hi < 0);
 141 }
 142
 143
 144 /**************************************************************************
 145                         m 6 4 M a c
 146 ** Mixed precision multiply and accumulate primitive for number building.
 147 ** Multiplies DPUNS u by FICL_UNS mul and adds FICL_UNS add. Mul is typically
 148 ** the numeric base, and add represents a digit to be appended to the
 149 ** growing number.
 150 ** Returns the result of the operation
 151 **************************************************************************/
 152 DPUNS m64Mac(DPUNS u, FICL_UNS mul, FICL_UNS add)
 153 {
 154     DPUNS resultLo = ficlLongMul(u.lo, mul);
 155     DPUNS resultHi = ficlLongMul(u.hi, mul);
 156     resultLo.hi += resultHi.lo;
 157     resultHi.lo = resultLo.lo + add;
 158
 159     if (resultHi.lo < resultLo.lo)
 160         resultLo.hi++;
 161
 162     resultLo.lo = resultHi.lo;
 163
 164     return resultLo;
 165 }
 166
 167
 168 /**************************************************************************
 169                         m 6 4 M u l I
 170 ** Multiplies a pair of FICL_INTs and returns an DPINT result.
 171 **************************************************************************/
 172 DPINT m64MulI(FICL_INT x, FICL_INT y)
 173 {
 174     DPUNS prod;
 175     int sign = 1;
 176
 177     if (x < 0)
 178     {
 179         sign = -sign;
 180         x = -x;
 181     }
 182
 183     if (y < 0)
 184     {
 185         sign = -sign;
 186         y = -y;
 187     }
 188
 189     prod = ficlLongMul(x, y);
 190     if (sign > 0)
 191         return m64CastUI(prod);
 192     else
 193         return m64Negate(m64CastUI(prod));
 194 }
 195
 196
 197 /**************************************************************************
 198                         m 6 4 N e g a t e
 199 ** Negates an DPINT by complementing and incrementing.
 200 **************************************************************************/
 201 DPINT m64Negate(DPINT x)
 202 {
 203     x.hi = ~x.hi;
 204     x.lo = ~x.lo;
 205     x.lo ++;
 206     if (x.lo == 0)
 207         x.hi++;
 208
 209     return x;
 210 }
 211
 212
 213 /**************************************************************************
 214                         m 6 4 P u s h
 215 ** Push an DPINT onto the specified stack in the order required
 216 ** by ANS Forth (most significant cell on top)
 217 ** These should probably be macros...
 218 **************************************************************************/
 219 void  i64Push(FICL_STACK *pStack, DPINT i64)
 220 {
 221     stackPushINT(pStack, i64.lo);
 222     stackPushINT(pStack, i64.hi);
 223     return;
 224 }
 225
 226 void  u64Push(FICL_STACK *pStack, DPUNS u64)
 227 {
 228     stackPushINT(pStack, u64.lo);
 229     stackPushINT(pStack, u64.hi);
 230     return;
 231 }
 232
 233
 234 /**************************************************************************
 235                         m 6 4 P o p
 236 ** Pops an DPINT off the stack in the order required by ANS Forth
 237 ** (most significant cell on top)
 238 ** These should probably be macros...
 239 **************************************************************************/
 240 DPINT i64Pop(FICL_STACK *pStack)
 241 {
 242     DPINT ret;
 243     ret.hi = stackPopINT(pStack);
 244     ret.lo = stackPopINT(pStack);
 245     return ret;
 246 }
 247
 248 DPUNS u64Pop(FICL_STACK *pStack)
 249 {
 250     DPUNS ret;
 251     ret.hi = stackPopINT(pStack);
 252     ret.lo = stackPopINT(pStack);
 253     return ret;
 254 }
 255
 256
 257 /**************************************************************************
 258                         m 6 4 S y m m e t r i c D i v
 259 ** Divide an DPINT by a FICL_INT and return a FICL_INT quotient and a
 260 ** FICL_INT remainder. The absolute values of quotient and remainder are not
 261 ** affected by the signs of the numerator and denominator (the operation
 262 ** is symmetric on the number line)
 263 **************************************************************************/
 264 INTQR m64SymmetricDivI(DPINT num, FICL_INT den)
 265 {
 266     INTQR qr;
 267     UNSQR uqr;
 268     int signRem = 1;
 269     int signQuot = 1;
 270
 271     if (m64IsNegative(num))
 272     {
 273         num = m64Negate(num);
 274         signRem  = -signRem;
 275         signQuot = -signQuot;
 276     }
 277
 278     if (den < 0)
 279     {
 280         den      = -den;
 281         signQuot = -signQuot;
 282     }
 283
 284     uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
 285     qr = m64CastQRUI(uqr);
 286     if (signRem < 0)
 287         qr.rem = -qr.rem;
 288
 289     if (signQuot < 0)
 290         qr.quot = -qr.quot;
 291
 292     return qr;
 293 }
 294
 295
 296 /**************************************************************************
 297                         m 6 4 U M o d
 298 ** Divides a DPUNS by base (an UNS16) and returns an UNS16 remainder.
 299 ** Writes the quotient back to the original DPUNS as a side effect.
 300 ** This operation is typically used to convert an DPUNS to a text string
 301 ** in any base. See words.c:numberSignS, for example.
 302 ** Mechanics: performs 4 ficlLongDivs, each of which produces 16 bits
 303 ** of the quotient. C does not provide a way to divide an FICL_UNS by an
 304 ** UNS16 and get an FICL_UNS quotient (ldiv is closest, but it's signed,
 305 ** unfortunately), so I've used ficlLongDiv.
 306 **************************************************************************/
 307 #if (BITS_PER_CELL == 32)
 308
 309 #define UMOD_SHIFT 16
 310 #define UMOD_MASK 0x0000ffff
 311
 312 #elif (BITS_PER_CELL == 64)
 313
 314 #define UMOD_SHIFT 32
 315 #define UMOD_MASK 0x00000000ffffffff
 316
 317 #endif
 318
 319 UNS16 m64UMod(DPUNS *pUD, UNS16 base)
 320 {
 321     DPUNS ud;
 322     UNSQR qr;
 323     DPUNS result;
 324
 325     result.hi = result.lo = 0;
 326
 327     ud.hi = 0;
 328     ud.lo = pUD->hi >> UMOD_SHIFT;
 329     qr = ficlLongDiv(ud, (FICL_UNS)base);
 330     result.hi = qr.quot << UMOD_SHIFT;
 331
 332     ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->hi & UMOD_MASK);
 333     qr = ficlLongDiv(ud, (FICL_UNS)base);
 334     result.hi |= qr.quot & UMOD_MASK;
 335
 336     ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo >> UMOD_SHIFT);
 337     qr = ficlLongDiv(ud, (FICL_UNS)base);
 338     result.lo = qr.quot << UMOD_SHIFT;
 339
 340     ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo & UMOD_MASK);
 341     qr = ficlLongDiv(ud, (FICL_UNS)base);
 342     result.lo |= qr.quot & UMOD_MASK;
 343
 344     *pUD = result;
 345
 346     return (UNS16)(qr.rem);
 347 }
 348
 349
 350 /**************************************************************************
 351 ** Contributed by
 352 ** Michael A. Gauland   gaulandm@mdhost.cse.tek.com
 353 **************************************************************************/
 354 #if PORTABLE_LONGMULDIV != 0
 355 /**************************************************************************
 356                         m 6 4 A d d
 357 **
 358 **************************************************************************/
 359 DPUNS m64Add(DPUNS x, DPUNS y)
 360 {
 361     DPUNS result;
 362     int carry;
 363
 364     result.hi = x.hi + y.hi;
 365     result.lo = x.lo + y.lo;
 366
 367
 368     carry  = ((x.lo | y.lo) & CELL_HI_BIT) && !(result.lo & CELL_HI_BIT);
 369     carry |= ((x.lo & y.lo) & CELL_HI_BIT);
 370
 371     if (carry)
 372     {
 373         result.hi++;
 374     }
 375
 376     return result;
 377 }
 378
 379
 380 /**************************************************************************
 381                         m 6 4 S u b
 382 **
 383 **************************************************************************/
 384 DPUNS m64Sub(DPUNS x, DPUNS y)
 385 {
 386     DPUNS result;
 387
 388     result.hi = x.hi - y.hi;
 389     result.lo = x.lo - y.lo;
 390
 391     if (x.lo < y.lo)
 392     {
 393         result.hi--;
 394     }
 395
 396     return result;
 397 }
 398
 399
 400 /**************************************************************************
 401                         m 6 4 A S L
 402 ** 64 bit left shift
 403 **************************************************************************/
 404 DPUNS m64ASL( DPUNS x )
 405 {
 406     DPUNS result;
 407
 408     result.hi = x.hi << 1;
 409     if (x.lo & CELL_HI_BIT)
 410     {
 411         result.hi++;
 412     }
 413
 414     result.lo = x.lo << 1;
 415
 416     return result;
 417 }
 418
 419
 420 /**************************************************************************
 421                         m 6 4 A S R
 422 ** 64 bit right shift (unsigned - no sign extend)
 423 **************************************************************************/
 424 DPUNS m64ASR( DPUNS x )
 425 {
 426     DPUNS result;
 427
 428     result.lo = x.lo >> 1;
 429     if (x.hi & 1)
 430     {
 431         result.lo |= CELL_HI_BIT;
 432     }
 433
 434     result.hi = x.hi >> 1;
 435     return result;
 436 }
 437
 438
 439 /**************************************************************************
 440                         m 6 4 O r
 441 ** 64 bit bitwise OR
 442 **************************************************************************/
 443 DPUNS m64Or( DPUNS x, DPUNS y )
 444 {
 445     DPUNS result;
 446
 447     result.hi = x.hi | y.hi;
 448     result.lo = x.lo | y.lo;
 449
 450     return result;
 451 }
 452
 453
 454 /**************************************************************************
 455                         m 6 4 C o m p a r e
 456 ** Return -1 if x < y; 0 if x==y, and 1 if x > y.
 457 **************************************************************************/
 458 int m64Compare(DPUNS x, DPUNS y)
 459 {
 460     int result;
 461
 462     if (x.hi > y.hi)
 463     {
 464         result = +1;
 465     }
 466     else if (x.hi < y.hi)
 467     {
 468         result = -1;
 469     }
 470     else
 471     {
 472         /* High parts are equal */
 473         if (x.lo > y.lo)
 474         {
 475             result = +1;
 476         }
 477         else if (x.lo < y.lo)
 478         {
 479             result = -1;
 480         }
 481         else
 482         {
 483             result = 0;
 484         }
 485     }
 486
 487     return result;
 488 }
 489
 490
 491 /**************************************************************************
 492                         f i c l L o n g M u l
 493 ** Portable versions of ficlLongMul and ficlLongDiv in C
 494 ** Contributed by:
 495 ** Michael A. Gauland   gaulandm@mdhost.cse.tek.com
 496 **************************************************************************/
 497 DPUNS ficlLongMul(FICL_UNS x, FICL_UNS y)
 498 {
 499     DPUNS result = { 0, 0 };
 500     DPUNS addend;
 501
 502     addend.lo = y;
 503     addend.hi = 0; /* No sign extension--arguments are unsigned */
 504
 505     while (x != 0)
 506     {
 507         if ( x & 1)
 508         {
 509             result = m64Add(result, addend);
 510         }
 511         x >>= 1;
 512         addend = m64ASL(addend);
 513     }
 514     return result;
 515 }
 516
 517
 518 /**************************************************************************
 519                         f i c l L o n g D i v
 520 ** Portable versions of ficlLongMul and ficlLongDiv in C
 521 ** Contributed by:
 522 ** Michael A. Gauland   gaulandm@mdhost.cse.tek.com
 523 **************************************************************************/
 524 UNSQR ficlLongDiv(DPUNS q, FICL_UNS y)
 525 {
 526     UNSQR result;
 527     DPUNS quotient;
 528     DPUNS subtrahend;
 529     DPUNS mask;
 530
 531     quotient.lo = 0;
 532     quotient.hi = 0;
 533
 534     subtrahend.lo = y;
 535     subtrahend.hi = 0;
 536
 537     mask.lo = 1;
 538     mask.hi = 0;
 539
 540     while ((m64Compare(subtrahend, q) < 0) &&
 541            (subtrahend.hi & CELL_HI_BIT) == 0)
 542     {
 543         mask = m64ASL(mask);
 544         subtrahend = m64ASL(subtrahend);
 545     }
 546
 547     while (mask.lo != 0 || mask.hi != 0)
 548     {
 549         if (m64Compare(subtrahend, q) <= 0)
 550         {
 551             q = m64Sub( q, subtrahend);
 552             quotient = m64Or(quotient, mask);
 553         }
 554         mask = m64ASR(mask);
 555         subtrahend = m64ASR(subtrahend);
 556     }
 557
 558     result.quot = quotient.lo;
 559     result.rem = q.lo;
 560     return result;
 561 }
 562
 563 #endif
 564