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 561