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**************************************************************************/ 55DPINT 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**************************************************************************/ 92INTQR 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**************************************************************************/ 135int 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**************************************************************************/ 149DPUNS 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**************************************************************************/ 169DPINT 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**************************************************************************/ 198DPINT 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**************************************************************************/ 216void i64Push(FICL_STACK *pStack, DPINT i64) 217{ 218 stackPushINT(pStack, i64.lo); 219 stackPushINT(pStack, i64.hi); 220 return; 221} 222 223void 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**************************************************************************/ 237DPINT i64Pop(FICL_STACK *pStack) 238{ 239 DPINT ret; 240 ret.hi = stackPopINT(pStack); 241 ret.lo = stackPopINT(pStack); 242 return ret; 243} 244 245DPUNS 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**************************************************************************/ 261INTQR 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 316UNS16 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**************************************************************************/ 356DPUNS 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**************************************************************************/ 381DPUNS 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**************************************************************************/ 401DPUNS 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**************************************************************************/ 421DPUNS 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**************************************************************************/ 440DPUNS 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**************************************************************************/ 455int 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**************************************************************************/ 494DPUNS 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**************************************************************************/ 521UNSQR 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