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 *******************************************************************/
   9
  10 /* $FreeBSD: src/sys/boot/ficl/math64.c,v 1.2 1999/09/29 04:43:06 dcs Exp $ */
  11 /* $DragonFly: src/sys/boot/ficl/math64.c,v 1.2 2003/06/17 04:28:17 dillon Exp $ */
  12
  13 #include "ficl.h"
  14 #include "math64.h"
  15
  16
  17 /**************************************************************************
  18                         m 6 4 A b s
  19 ** Returns the absolute value of an DPINT
  20 **************************************************************************/
  21 DPINT m64Abs(DPINT x)
  22 {
  23     if (m64IsNegative(x))
  24         x = m64Negate(x);
  25
  26     return x;
  27 }
  28
  29
  30 /**************************************************************************
  31                         m 6 4 F l o o r e d D i v I
  32 **
  33 ** FROM THE FORTH ANS...
  34 ** Floored division is integer division in which the remainder carries
  35 ** the sign of the divisor or is zero, and the quotient is rounded to
  36 ** its arithmetic floor. Symmetric division is integer division in which
  37 ** the remainder carries the sign of the dividend or is zero and the
  38 ** quotient is the mathematical quotient rounded towards zero or
  39 ** truncated. Examples of each are shown in tables 3.3 and 3.4.
  40 **
  41 ** Table 3.3 - Floored Division Example
  42 ** Dividend        Divisor Remainder       Quotient
  43 ** --------        ------- ---------       --------
  44 **  10                7       3                1
  45 ** -10                7       4               -2
  46 **  10               -7      -4               -2
  47 ** -10               -7      -3                1
  48 **
  49 **
  50 ** Table 3.4 - Symmetric Division Example
  51 ** Dividend        Divisor Remainder       Quotient
  52 ** --------        ------- ---------       --------
  53 **  10                7       3                1
  54 ** -10                7      -3               -1
  55 **  10               -7       3               -1
  56 ** -10               -7      -3                1
  57 **************************************************************************/
  58 INTQR m64FlooredDivI(DPINT num, FICL_INT den)
  59 {
  60     INTQR qr;
  61     UNSQR uqr;
  62     int signRem = 1;
  63     int signQuot = 1;
  64
  65     if (m64IsNegative(num))
  66     {
  67         num = m64Negate(num);
  68         signQuot = -signQuot;
  69     }
  70
  71     if (den < 0)
  72     {
  73         den      = -den;
  74         signRem  = -signRem;
  75         signQuot = -signQuot;
  76     }
  77
  78     uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
  79     qr = m64CastQRUI(uqr);
  80     if (signQuot < 0)
  81     {
  82         qr.quot = -qr.quot;
  83         if (qr.rem != 0)
  84         {
  85             qr.quot--;
  86             qr.rem = den - qr.rem;
  87         }
  88     }
  89
  90     if (signRem < 0)
  91         qr.rem = -qr.rem;
  92
  93     return qr;
  94 }
  95
  96
  97 /**************************************************************************
  98                         m 6 4 I s N e g a t i v e
  99 ** Returns TRUE if the specified DPINT has its sign bit set.
 100 **************************************************************************/
 101 int m64IsNegative(DPINT x)
 102 {
 103     return (x.hi < 0);
 104 }
 105
 106
 107 /**************************************************************************
 108                         m 6 4 M a c
 109 ** Mixed precision multiply and accumulate primitive for number building.
 110 ** Multiplies DPUNS u by FICL_UNS mul and adds FICL_UNS add. Mul is typically
 111 ** the numeric base, and add represents a digit to be appended to the
 112 ** growing number.
 113 ** Returns the result of the operation
 114 **************************************************************************/
 115 DPUNS m64Mac(DPUNS u, FICL_UNS mul, FICL_UNS add)
 116 {
 117     DPUNS resultLo = ficlLongMul(u.lo, mul);
 118     DPUNS resultHi = ficlLongMul(u.hi, mul);
 119     resultLo.hi += resultHi.lo;
 120     resultHi.lo = resultLo.lo + add;
 121
 122     if (resultHi.lo < resultLo.lo)
 123         resultLo.hi++;
 124
 125     resultLo.lo = resultHi.lo;
 126
 127     return resultLo;
 128 }
 129
 130
 131 /**************************************************************************
 132                         m 6 4 M u l I
 133 ** Multiplies a pair of FICL_INTs and returns an DPINT result.
 134 **************************************************************************/
 135 DPINT m64MulI(FICL_INT x, FICL_INT y)
 136 {
 137     DPUNS prod;
 138     int sign = 1;
 139
 140     if (x < 0)
 141     {
 142         sign = -sign;
 143         x = -x;
 144     }
 145
 146     if (y < 0)
 147     {
 148         sign = -sign;
 149         y = -y;
 150     }
 151
 152     prod = ficlLongMul(x, y);
 153     if (sign > 0)
 154         return m64CastUI(prod);
 155     else
 156         return m64Negate(m64CastUI(prod));
 157 }
 158
 159
 160 /**************************************************************************
 161                         m 6 4 N e g a t e
 162 ** Negates an DPINT by complementing and incrementing.
 163 **************************************************************************/
 164 DPINT m64Negate(DPINT x)
 165 {
 166     x.hi = ~x.hi;
 167     x.lo = ~x.lo;
 168     x.lo ++;
 169     if (x.lo == 0)
 170         x.hi++;
 171
 172     return x;
 173 }
 174
 175
 176 /**************************************************************************
 177                         m 6 4 P u s h
 178 ** Push an DPINT onto the specified stack in the order required
 179 ** by ANS Forth (most significant cell on top)
 180 ** These should probably be macros...
 181 **************************************************************************/
 182 void  i64Push(FICL_STACK *pStack, DPINT i64)
 183 {
 184     stackPushINT(pStack, i64.lo);
 185     stackPushINT(pStack, i64.hi);
 186     return;
 187 }
 188
 189 void  u64Push(FICL_STACK *pStack, DPUNS u64)
 190 {
 191     stackPushINT(pStack, u64.lo);
 192     stackPushINT(pStack, u64.hi);
 193     return;
 194 }
 195
 196
 197 /**************************************************************************
 198                         m 6 4 P o p
 199 ** Pops an DPINT off the stack in the order required by ANS Forth
 200 ** (most significant cell on top)
 201 ** These should probably be macros...
 202 **************************************************************************/
 203 DPINT i64Pop(FICL_STACK *pStack)
 204 {
 205     DPINT ret;
 206     ret.hi = stackPopINT(pStack);
 207     ret.lo = stackPopINT(pStack);
 208     return ret;
 209 }
 210
 211 DPUNS u64Pop(FICL_STACK *pStack)
 212 {
 213     DPUNS ret;
 214     ret.hi = stackPopINT(pStack);
 215     ret.lo = stackPopINT(pStack);
 216     return ret;
 217 }
 218
 219
 220 /**************************************************************************
 221                         m 6 4 S y m m e t r i c D i v
 222 ** Divide an DPINT by a FICL_INT and return a FICL_INT quotient and a
 223 ** FICL_INT remainder. The absolute values of quotient and remainder are not
 224 ** affected by the signs of the numerator and denominator (the operation
 225 ** is symmetric on the number line)
 226 **************************************************************************/
 227 INTQR m64SymmetricDivI(DPINT num, FICL_INT den)
 228 {
 229     INTQR qr;
 230     UNSQR uqr;
 231     int signRem = 1;
 232     int signQuot = 1;
 233
 234     if (m64IsNegative(num))
 235     {
 236         num = m64Negate(num);
 237         signRem  = -signRem;
 238         signQuot = -signQuot;
 239     }
 240
 241     if (den < 0)
 242     {
 243         den      = -den;
 244         signQuot = -signQuot;
 245     }
 246
 247     uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
 248     qr = m64CastQRUI(uqr);
 249     if (signRem < 0)
 250         qr.rem = -qr.rem;
 251
 252     if (signQuot < 0)
 253         qr.quot = -qr.quot;
 254
 255     return qr;
 256 }
 257
 258
 259 /**************************************************************************
 260                         m 6 4 U M o d
 261 ** Divides a DPUNS by base (an UNS16) and returns an UNS16 remainder.
 262 ** Writes the quotient back to the original DPUNS as a side effect.
 263 ** This operation is typically used to convert an DPUNS to a text string
 264 ** in any base. See words.c:numberSignS, for example.
 265 ** Mechanics: performs 4 ficlLongDivs, each of which produces 16 bits
 266 ** of the quotient. C does not provide a way to divide an FICL_UNS by an
 267 ** UNS16 and get an FICL_UNS quotient (ldiv is closest, but it's signed,
 268 ** unfortunately), so I've used ficlLongDiv.
 269 **************************************************************************/
 270 #if (BITS_PER_CELL == 32)
 271
 272 #define UMOD_SHIFT 16
 273 #define UMOD_MASK 0x0000ffff
 274
 275 #elif (BITS_PER_CELL == 64)
 276
 277 #define UMOD_SHIFT 32
 278 #define UMOD_MASK 0x00000000ffffffff
 279
 280 #endif
 281
 282 UNS16 m64UMod(DPUNS *pUD, UNS16 base)
 283 {
 284     DPUNS ud;
 285     UNSQR qr;
 286     DPUNS result;
 287
 288     result.hi = result.lo = 0;
 289
 290     ud.hi = 0;
 291     ud.lo = pUD->hi >> UMOD_SHIFT;
 292     qr = ficlLongDiv(ud, (FICL_UNS)base);
 293     result.hi = qr.quot << UMOD_SHIFT;
 294
 295     ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->hi & UMOD_MASK);
 296     qr = ficlLongDiv(ud, (FICL_UNS)base);
 297     result.hi |= qr.quot & UMOD_MASK;
 298
 299     ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo >> UMOD_SHIFT);
 300     qr = ficlLongDiv(ud, (FICL_UNS)base);
 301     result.lo = qr.quot << UMOD_SHIFT;
 302
 303     ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo & UMOD_MASK);
 304     qr = ficlLongDiv(ud, (FICL_UNS)base);
 305     result.lo |= qr.quot & UMOD_MASK;
 306
 307     *pUD = result;
 308
 309     return (UNS16)(qr.rem);
 310 }
 311
 312
 313 /**************************************************************************
 314 ** Contributed by
 315 ** Michael A. Gauland   gaulandm@mdhost.cse.tek.com
 316 **************************************************************************/
 317 #if PORTABLE_LONGMULDIV != 0
 318 /**************************************************************************
 319                         m 6 4 A d d
 320 **
 321 **************************************************************************/
 322 DPUNS m64Add(DPUNS x, DPUNS y)
 323 {
 324     DPUNS result;
 325     int carry;
 326
 327     result.hi = x.hi + y.hi;
 328     result.lo = x.lo + y.lo;
 329
 330
 331     carry  = ((x.lo | y.lo) & CELL_HI_BIT) && !(result.lo & CELL_HI_BIT);
 332     carry |= ((x.lo & y.lo) & CELL_HI_BIT);
 333
 334     if (carry)
 335     {
 336         result.hi++;
 337     }
 338
 339     return result;
 340 }
 341
 342
 343 /**************************************************************************
 344                         m 6 4 S u b
 345 **
 346 **************************************************************************/
 347 DPUNS m64Sub(DPUNS x, DPUNS y)
 348 {
 349     DPUNS result;
 350
 351     result.hi = x.hi - y.hi;
 352     result.lo = x.lo - y.lo;
 353
 354     if (x.lo < y.lo)
 355     {
 356         result.hi--;
 357     }
 358
 359     return result;
 360 }
 361
 362
 363 /**************************************************************************
 364                         m 6 4 A S L
 365 ** 64 bit left shift
 366 **************************************************************************/
 367 DPUNS m64ASL( DPUNS x )
 368 {
 369     DPUNS result;
 370
 371     result.hi = x.hi << 1;
 372     if (x.lo & CELL_HI_BIT)
 373     {
 374         result.hi++;
 375     }
 376
 377     result.lo = x.lo << 1;
 378
 379     return result;
 380 }
 381
 382
 383 /**************************************************************************
 384                         m 6 4 A S R
 385 ** 64 bit right shift (unsigned - no sign extend)
 386 **************************************************************************/
 387 DPUNS m64ASR( DPUNS x )
 388 {
 389     DPUNS result;
 390
 391     result.lo = x.lo >> 1;
 392     if (x.hi & 1)
 393     {
 394         result.lo |= CELL_HI_BIT;
 395     }
 396
 397     result.hi = x.hi >> 1;
 398     return result;
 399 }
 400
 401
 402 /**************************************************************************
 403                         m 6 4 O r
 404 ** 64 bit bitwise OR
 405 **************************************************************************/
 406 DPUNS m64Or( DPUNS x, DPUNS y )
 407 {
 408     DPUNS result;
 409
 410     result.hi = x.hi | y.hi;
 411     result.lo = x.lo | y.lo;
 412
 413     return result;
 414 }
 415
 416
 417 /**************************************************************************
 418                         m 6 4 C o m p a r e
 419 ** Return -1 if x < y; 0 if x==y, and 1 if x > y.
 420 **************************************************************************/
 421 int m64Compare(DPUNS x, DPUNS y)
 422 {
 423     int result;
 424
 425     if (x.hi > y.hi)
 426     {
 427         result = +1;
 428     }
 429     else if (x.hi < y.hi)
 430     {
 431         result = -1;
 432     }
 433     else
 434     {
 435         /* High parts are equal */
 436         if (x.lo > y.lo)
 437         {
 438             result = +1;
 439         }
 440         else if (x.lo < y.lo)
 441         {
 442             result = -1;
 443         }
 444         else
 445         {
 446             result = 0;
 447         }
 448     }
 449
 450     return result;
 451 }
 452
 453
 454 /**************************************************************************
 455                         f i c l L o n g M u l
 456 ** Portable versions of ficlLongMul and ficlLongDiv in C
 457 ** Contributed by:
 458 ** Michael A. Gauland   gaulandm@mdhost.cse.tek.com
 459 **************************************************************************/
 460 DPUNS ficlLongMul(FICL_UNS x, FICL_UNS y)
 461 {
 462     DPUNS result = { 0, 0 };
 463     DPUNS addend;
 464
 465     addend.lo = y;
 466     addend.hi = 0; /* No sign extension--arguments are unsigned */
 467
 468     while (x != 0)
 469     {
 470         if ( x & 1)
 471         {
 472             result = m64Add(result, addend);
 473         }
 474         x >>= 1;
 475         addend = m64ASL(addend);
 476     }
 477     return result;
 478 }
 479
 480
 481 /**************************************************************************
 482                         f i c l L o n g D i v
 483 ** Portable versions of ficlLongMul and ficlLongDiv in C
 484 ** Contributed by:
 485 ** Michael A. Gauland   gaulandm@mdhost.cse.tek.com
 486 **************************************************************************/
 487 UNSQR ficlLongDiv(DPUNS q, FICL_UNS y)
 488 {
 489     UNSQR result;
 490     DPUNS quotient;
 491     DPUNS subtrahend;
 492     DPUNS mask;
 493
 494     quotient.lo = 0;
 495     quotient.hi = 0;
 496
 497     subtrahend.lo = y;
 498     subtrahend.hi = 0;
 499
 500     mask.lo = 1;
 501     mask.hi = 0;
 502
 503     while ((m64Compare(subtrahend, q) < 0) &&
 504            (subtrahend.hi & CELL_HI_BIT) == 0)
 505     {
 506         mask = m64ASL(mask);
 507         subtrahend = m64ASL(subtrahend);
 508     }
 509
 510     while (mask.lo != 0 || mask.hi != 0)
 511     {
 512         if (m64Compare(subtrahend, q) <= 0)
 513         {
 514             q = m64Sub( q, subtrahend);
 515             quotient = m64Or(quotient, mask);
 516         }
 517         mask = m64ASR(mask);
 518         subtrahend = m64ASR(subtrahend);
 519     }
 520
 521     result.quot = quotient.lo;
 522     result.rem = q.lo;
 523     return result;
 524 }
 525
 526 #endif
 527