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