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**************************************************************************/ 54DPINT 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**************************************************************************/ 91INTQR 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**************************************************************************/ 134int 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**************************************************************************/ 148DPUNS 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**************************************************************************/ 168DPINT 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**************************************************************************/ 197DPINT 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**************************************************************************/ 215void i64Push(FICL_STACK *pStack, DPINT i64) 216{ 217 stackPushINT(pStack, i64.lo); 218 stackPushINT(pStack, i64.hi); 219 return; 220} 221 222void 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**************************************************************************/ 236DPINT i64Pop(FICL_STACK *pStack) 237{ 238 DPINT ret; 239 ret.hi = stackPopINT(pStack); 240 ret.lo = stackPopINT(pStack); 241 return ret; 242} 243 244DPUNS 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**************************************************************************/ 260INTQR 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 315UNS16 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**************************************************************************/ 355DPUNS 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**************************************************************************/ 380DPUNS 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**************************************************************************/ 400DPUNS 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**************************************************************************/ 420DPUNS 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**************************************************************************/ 439DPUNS 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**************************************************************************/ 454int 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**************************************************************************/ 493DPUNS 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**************************************************************************/ 520UNSQR 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