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