1 /* 2 * m a t h 6 4 . c 3 * Forth Inspired Command Language - 64 bit math support routines 4 * Authors: Michael A. Gauland (gaulandm@mdhost.cse.tek.com) 5 * Larry Hastings (larry@hastings.org) 6 * John Sadler (john_sadler@alum.mit.edu) 7 * Created: 25 January 1998 8 * Rev 2.03: Support for 128 bit DP math. This file really ouught to 9 * be renamed! 10 * $Id: double.c,v 1.2 2010/09/12 15:18:07 asau Exp $ 11 */ 12 /* 13 * Copyright (c) 1997-2001 John Sadler (john_sadler@alum.mit.edu) 14 * All rights reserved. 15 * 16 * Get the latest Ficl release at http://ficl.sourceforge.net 17 * 18 * I am interested in hearing from anyone who uses Ficl. If you have 19 * a problem, a success story, a defect, an enhancement request, or 20 * if you would like to contribute to the Ficl release, please 21 * contact me by email at the address above. 22 * 23 * L I C E N S E and D I S C L A I M E R 24 * 25 * Redistribution and use in source and binary forms, with or without 26 * modification, are permitted provided that the following conditions 27 * are met: 28 * 1. Redistributions of source code must retain the above copyright 29 * notice, this list of conditions and the following disclaimer. 30 * 2. Redistributions in binary form must reproduce the above copyright 31 * notice, this list of conditions and the following disclaimer in the 32 * documentation and/or other materials provided with the distribution. 33 * 34 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 35 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 36 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 37 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 38 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 39 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 40 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 41 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 42 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 43 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 44 * SUCH DAMAGE. 45 */ 46 47 #include "ficl.h" 48 49 #if FICL_PLATFORM_HAS_2INTEGER 50 ficl2UnsignedQR 51 ficl2UnsignedDivide(ficl2Unsigned q, ficlUnsigned y) 52 { 53 ficl2UnsignedQR result; 54 55 result.quotient = q / y; 56 /* 57 * Once we have the quotient, it's cheaper to calculate the 58 * remainder this way than with % (mod). --lch 59 */ 60 result.remainder = (ficlInteger)(q - (result.quotient * y)); 61 62 return (result); 63 } 64 65 #else /* FICL_PLATFORM_HAS_2INTEGER */ 66 67 #define FICL_CELL_HIGH_BIT ((uintmax_t)1 << (FICL_BITS_PER_CELL-1)) 68 #define UMOD_SHIFT (FICL_BITS_PER_CELL / 2) 69 #define UMOD_MASK ((1L << (FICL_BITS_PER_CELL / 2)) - 1) 70 71 /* 72 * ficl2IntegerIsNegative 73 * Returns TRUE if the specified ficl2Unsigned has its sign bit set. 74 */ 75 int 76 ficl2IntegerIsNegative(ficl2Integer x) 77 { 78 return (x.high < 0); 79 } 80 81 /* 82 * ficl2IntegerNegate 83 * Negates an ficl2Unsigned by complementing and incrementing. 84 */ 85 ficl2Integer 86 ficl2IntegerNegate(ficl2Integer x) 87 { 88 x.high = ~x.high; 89 x.low = ~x.low; 90 x.low ++; 91 if (x.low == 0) 92 x.high++; 93 94 return (x); 95 } 96 97 /* 98 * ficl2UnsignedMultiplyAccumulate 99 * Mixed precision multiply and accumulate primitive for number building. 100 * Multiplies ficl2Unsigned u by ficlUnsigned mul and adds ficlUnsigned add. 101 * Mul is typically the numeric base, and add represents a digit to be 102 * appended to the growing number. 103 * Returns the result of the operation 104 */ 105 ficl2Unsigned 106 ficl2UnsignedMultiplyAccumulate(ficl2Unsigned u, ficlUnsigned mul, 107 ficlUnsigned add) 108 { 109 ficl2Unsigned resultLo = ficl2UnsignedMultiply(u.low, mul); 110 ficl2Unsigned resultHi = ficl2UnsignedMultiply(u.high, mul); 111 resultLo.high += resultHi.low; 112 resultHi.low = resultLo.low + add; 113 114 if (resultHi.low < resultLo.low) 115 resultLo.high++; 116 117 resultLo.low = resultHi.low; 118 119 return (resultLo); 120 } 121 122 /* 123 * ficl2IntegerMultiply 124 * Multiplies a pair of ficlIntegers and returns an ficl2Integer result. 125 */ 126 ficl2Integer 127 ficl2IntegerMultiply(ficlInteger x, ficlInteger y) 128 { 129 ficl2Unsigned prod; 130 ficl2Integer result; 131 int sign = 1; 132 133 if (x < 0) { 134 sign = -sign; 135 x = -x; 136 } 137 138 if (y < 0) { 139 sign = -sign; 140 y = -y; 141 } 142 143 prod = ficl2UnsignedMultiply(x, y); 144 FICL_2INTEGER_SET(FICL_2UNSIGNED_GET_HIGH(prod), 145 FICL_2UNSIGNED_GET_LOW(prod), result); 146 if (sign > 0) 147 return (result); 148 else 149 return (ficl2IntegerNegate(result)); 150 } 151 152 ficl2Integer 153 ficl2IntegerDecrement(ficl2Integer x) 154 { 155 if (x.low == INTMAX_MIN) 156 x.high--; 157 x.low--; 158 159 return (x); 160 } 161 162 ficl2Unsigned 163 ficl2UnsignedAdd(ficl2Unsigned x, ficl2Unsigned y) 164 { 165 ficl2Unsigned result; 166 167 result.high = x.high + y.high; 168 result.low = x.low + y.low; 169 170 if (result.low < y.low) 171 result.high++; 172 173 return (result); 174 } 175 176 /* 177 * ficl2UnsignedMultiply 178 * Contributed by: 179 * Michael A. Gauland gaulandm@mdhost.cse.tek.com 180 */ 181 ficl2Unsigned 182 ficl2UnsignedMultiply(ficlUnsigned x, ficlUnsigned y) 183 { 184 ficl2Unsigned result = { 0, 0 }; 185 ficl2Unsigned addend; 186 187 addend.low = y; 188 addend.high = 0; /* No sign extension--arguments are unsigned */ 189 190 while (x != 0) { 191 if (x & 1) { 192 result = ficl2UnsignedAdd(result, addend); 193 } 194 x >>= 1; 195 addend = ficl2UnsignedArithmeticShiftLeft(addend); 196 } 197 return (result); 198 } 199 200 /* 201 * ficl2UnsignedSubtract 202 */ 203 ficl2Unsigned 204 ficl2UnsignedSubtract(ficl2Unsigned x, ficl2Unsigned y) 205 { 206 ficl2Unsigned result; 207 208 result.high = x.high - y.high; 209 result.low = x.low - y.low; 210 211 if (x.low < y.low) { 212 result.high--; 213 } 214 215 return (result); 216 } 217 218 /* 219 * ficl2UnsignedArithmeticShiftLeft 220 * 64 bit left shift 221 */ 222 ficl2Unsigned 223 ficl2UnsignedArithmeticShiftLeft(ficl2Unsigned x) 224 { 225 ficl2Unsigned result; 226 227 result.high = x.high << 1; 228 if (x.low & FICL_CELL_HIGH_BIT) { 229 result.high++; 230 } 231 232 result.low = x.low << 1; 233 234 return (result); 235 } 236 237 /* 238 * ficl2UnsignedArithmeticShiftRight 239 * 64 bit right shift (unsigned - no sign extend) 240 */ 241 ficl2Unsigned 242 ficl2UnsignedArithmeticShiftRight(ficl2Unsigned x) 243 { 244 ficl2Unsigned result; 245 246 result.low = x.low >> 1; 247 if (x.high & 1) { 248 result.low |= FICL_CELL_HIGH_BIT; 249 } 250 251 result.high = x.high >> 1; 252 return (result); 253 } 254 255 /* 256 * ficl2UnsignedOr 257 * 64 bit bitwise OR 258 */ 259 ficl2Unsigned 260 ficl2UnsignedOr(ficl2Unsigned x, ficl2Unsigned y) 261 { 262 ficl2Unsigned result; 263 264 result.high = x.high | y.high; 265 result.low = x.low | y.low; 266 267 return (result); 268 } 269 270 /* 271 * ficl2UnsignedCompare 272 * Return -1 if x < y; 0 if x==y, and 1 if x > y. 273 */ 274 int 275 ficl2UnsignedCompare(ficl2Unsigned x, ficl2Unsigned y) 276 { 277 if (x.high > y.high) 278 return (1); 279 if (x.high < y.high) 280 return (-1); 281 282 /* High parts are equal */ 283 284 if (x.low > y.low) 285 return (1); 286 else if (x.low < y.low) 287 return (-1); 288 289 return (0); 290 } 291 292 /* 293 * ficl2UnsignedDivide 294 * Portable versions of ficl2Multiply and ficl2Divide in C 295 * Contributed by: 296 * Michael A. Gauland gaulandm@mdhost.cse.tek.com 297 */ 298 ficl2UnsignedQR 299 ficl2UnsignedDivide(ficl2Unsigned q, ficlUnsigned y) 300 { 301 ficl2UnsignedQR result; 302 ficl2Unsigned quotient; 303 ficl2Unsigned subtrahend; 304 ficl2Unsigned mask; 305 306 quotient.low = 0; 307 quotient.high = 0; 308 309 subtrahend.low = y; 310 subtrahend.high = 0; 311 312 mask.low = 1; 313 mask.high = 0; 314 315 while ((ficl2UnsignedCompare(subtrahend, q) < 0) && 316 (subtrahend.high & FICL_CELL_HIGH_BIT) == 0) { 317 mask = ficl2UnsignedArithmeticShiftLeft(mask); 318 subtrahend = ficl2UnsignedArithmeticShiftLeft(subtrahend); 319 } 320 321 while (mask.low != 0 || mask.high != 0) { 322 if (ficl2UnsignedCompare(subtrahend, q) <= 0) { 323 q = ficl2UnsignedSubtract(q, subtrahend); 324 quotient = ficl2UnsignedOr(quotient, mask); 325 } 326 mask = ficl2UnsignedArithmeticShiftRight(mask); 327 subtrahend = ficl2UnsignedArithmeticShiftRight(subtrahend); 328 } 329 330 result.quotient = quotient; 331 result.remainder = q.low; 332 return (result); 333 } 334 #endif /* !FICL_PLATFORM_HAS_2INTEGER */ 335 336 /* 337 * ficl2IntegerDivideFloored 338 * 339 * FROM THE FORTH ANS... 340 * Floored division is integer division in which the remainder carries 341 * the sign of the divisor or is zero, and the quotient is rounded to 342 * its arithmetic floor. Symmetric division is integer division in which 343 * the remainder carries the sign of the dividend or is zero and the 344 * quotient is the mathematical quotient rounded towards zero or 345 * truncated. Examples of each are shown in tables 3.3 and 3.4. 346 * 347 * Table 3.3 - Floored Division Example 348 * Dividend Divisor Remainder Quotient 349 * -------- ------- --------- -------- 350 * 10 7 3 1 351 * -10 7 4 -2 352 * 10 -7 -4 -2 353 * -10 -7 -3 1 354 * 355 * 356 * Table 3.4 - Symmetric Division Example 357 * Dividend Divisor Remainder Quotient 358 * -------- ------- --------- -------- 359 * 10 7 3 1 360 * -10 7 -3 -1 361 * 10 -7 3 -1 362 * -10 -7 -3 1 363 */ 364 ficl2IntegerQR 365 ficl2IntegerDivideFloored(ficl2Integer num, ficlInteger den) 366 { 367 ficl2IntegerQR qr; 368 ficl2UnsignedQR uqr; 369 ficl2Unsigned u; 370 int signRem = 1; 371 int signQuot = 1; 372 373 if (ficl2IntegerIsNegative(num)) { 374 num = ficl2IntegerNegate(num); 375 signQuot = -signQuot; 376 } 377 378 if (den < 0) { 379 den = -den; 380 signRem = -signRem; 381 signQuot = -signQuot; 382 } 383 384 FICL_2UNSIGNED_SET(FICL_2UNSIGNED_GET_HIGH(num), 385 FICL_2UNSIGNED_GET_LOW(num), u); 386 uqr = ficl2UnsignedDivide(u, (ficlUnsigned)den); 387 qr = FICL_2UNSIGNEDQR_TO_2INTEGERQR(uqr); 388 if (signQuot < 0) { 389 qr.quotient = ficl2IntegerNegate(qr.quotient); 390 if (qr.remainder != 0) { 391 qr.quotient = ficl2IntegerDecrement(qr.quotient); 392 qr.remainder = den - qr.remainder; 393 } 394 } 395 396 if (signRem < 0) 397 qr.remainder = -qr.remainder; 398 399 return (qr); 400 } 401 402 /* 403 * ficl2IntegerDivideSymmetric 404 * Divide an ficl2Unsigned by a ficlInteger and return a ficlInteger quotient 405 * and a ficlInteger remainder. The absolute values of quotient and remainder 406 * are not affected by the signs of the numerator and denominator 407 * (the operation is symmetric on the number line) 408 */ 409 ficl2IntegerQR 410 ficl2IntegerDivideSymmetric(ficl2Integer num, ficlInteger den) 411 { 412 ficl2IntegerQR qr; 413 ficl2UnsignedQR uqr; 414 ficl2Unsigned u; 415 int signRem = 1; 416 int signQuot = 1; 417 418 if (ficl2IntegerIsNegative(num)) { 419 num = ficl2IntegerNegate(num); 420 signRem = -signRem; 421 signQuot = -signQuot; 422 } 423 424 if (den < 0) { 425 den = -den; 426 signQuot = -signQuot; 427 } 428 429 FICL_2UNSIGNED_SET(FICL_2UNSIGNED_GET_HIGH(num), 430 FICL_2UNSIGNED_GET_LOW(num), u); 431 uqr = ficl2UnsignedDivide(u, (ficlUnsigned)den); 432 qr = FICL_2UNSIGNEDQR_TO_2INTEGERQR(uqr); 433 if (signRem < 0) 434 qr.remainder = -qr.remainder; 435 436 if (signQuot < 0) 437 qr.quotient = ficl2IntegerNegate(qr.quotient); 438 439 return (qr); 440 } 441