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