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 562