1/* Common base code for the decNumber C Library. 2 Copyright (C) 2007-2015 Free Software Foundation, Inc. 3 Contributed by IBM Corporation. Author Mike Cowlishaw. 4 5 This file is part of GCC. 6 7 GCC is free software; you can redistribute it and/or modify it under 8 the terms of the GNU General Public License as published by the Free 9 Software Foundation; either version 3, or (at your option) any later 10 version. 11 12 GCC is distributed in the hope that it will be useful, but WITHOUT ANY 13 WARRANTY; without even the implied warranty of MERCHANTABILITY or 14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 15 for more details. 16 17Under Section 7 of GPL version 3, you are granted additional 18permissions described in the GCC Runtime Library Exception, version 193.1, as published by the Free Software Foundation. 20 21You should have received a copy of the GNU General Public License and 22a copy of the GCC Runtime Library Exception along with this program; 23see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 24<http://www.gnu.org/licenses/>. */ 25 26/* ------------------------------------------------------------------ */ 27/* decBasic.c -- common base code for Basic decimal types */ 28/* ------------------------------------------------------------------ */ 29/* This module comprises code that is shared between decDouble and */ 30/* decQuad (but not decSingle). The main arithmetic operations are */ 31/* here (Add, Subtract, Multiply, FMA, and Division operators). */ 32/* */ 33/* Unlike decNumber, parameterization takes place at compile time */ 34/* rather than at runtime. The parameters are set in the decDouble.c */ 35/* (etc.) files, which then include this one to produce the compiled */ 36/* code. The functions here, therefore, are code shared between */ 37/* multiple formats. */ 38/* */ 39/* This must be included after decCommon.c. */ 40/* ------------------------------------------------------------------ */ 41/* Names here refer to decFloat rather than to decDouble, etc., and */ 42/* the functions are in strict alphabetical order. */ 43 44/* The compile-time flags SINGLE, DOUBLE, and QUAD are set up in */ 45/* decCommon.c */ 46#if !defined(QUAD) 47 #error decBasic.c must be included after decCommon.c 48#endif 49#if SINGLE 50 #error Routines in decBasic.c are for decDouble and decQuad only 51#endif 52 53/* Private constants */ 54#define DIVIDE 0x80000000 /* Divide operations [as flags] */ 55#define REMAINDER 0x40000000 /* .. */ 56#define DIVIDEINT 0x20000000 /* .. */ 57#define REMNEAR 0x10000000 /* .. */ 58 59/* Private functions (local, used only by routines in this module) */ 60static decFloat *decDivide(decFloat *, const decFloat *, 61 const decFloat *, decContext *, uInt); 62static decFloat *decCanonical(decFloat *, const decFloat *); 63static void decFiniteMultiply(bcdnum *, uByte *, const decFloat *, 64 const decFloat *); 65static decFloat *decInfinity(decFloat *, const decFloat *); 66static decFloat *decInvalid(decFloat *, decContext *); 67static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *, 68 decContext *); 69static Int decNumCompare(const decFloat *, const decFloat *, Flag); 70static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *, 71 enum rounding, Flag); 72static uInt decToInt32(const decFloat *, decContext *, enum rounding, 73 Flag, Flag); 74 75/* ------------------------------------------------------------------ */ 76/* decCanonical -- copy a decFloat, making canonical */ 77/* */ 78/* result gets the canonicalized df */ 79/* df is the decFloat to copy and make canonical */ 80/* returns result */ 81/* */ 82/* This is exposed via decFloatCanonical for Double and Quad only. */ 83/* This works on specials, too; no error or exception is possible. */ 84/* ------------------------------------------------------------------ */ 85static decFloat * decCanonical(decFloat *result, const decFloat *df) { 86 uInt encode, precode, dpd; /* work */ 87 uInt inword, uoff, canon; /* .. */ 88 Int n; /* counter (down) */ 89 if (df!=result) *result=*df; /* effect copy if needed */ 90 if (DFISSPECIAL(result)) { 91 if (DFISINF(result)) return decInfinity(result, df); /* clean Infinity */ 92 /* is a NaN */ 93 DFWORD(result, 0)&=~ECONNANMASK; /* clear ECON except selector */ 94 if (DFISCCZERO(df)) return result; /* coefficient continuation is 0 */ 95 /* drop through to check payload */ 96 } 97 /* return quickly if the coefficient continuation is canonical */ 98 { /* declare block */ 99 #if DOUBLE 100 uInt sourhi=DFWORD(df, 0); 101 uInt sourlo=DFWORD(df, 1); 102 if (CANONDPDOFF(sourhi, 8) 103 && CANONDPDTWO(sourhi, sourlo, 30) 104 && CANONDPDOFF(sourlo, 20) 105 && CANONDPDOFF(sourlo, 10) 106 && CANONDPDOFF(sourlo, 0)) return result; 107 #elif QUAD 108 uInt sourhi=DFWORD(df, 0); 109 uInt sourmh=DFWORD(df, 1); 110 uInt sourml=DFWORD(df, 2); 111 uInt sourlo=DFWORD(df, 3); 112 if (CANONDPDOFF(sourhi, 4) 113 && CANONDPDTWO(sourhi, sourmh, 26) 114 && CANONDPDOFF(sourmh, 16) 115 && CANONDPDOFF(sourmh, 6) 116 && CANONDPDTWO(sourmh, sourml, 28) 117 && CANONDPDOFF(sourml, 18) 118 && CANONDPDOFF(sourml, 8) 119 && CANONDPDTWO(sourml, sourlo, 30) 120 && CANONDPDOFF(sourlo, 20) 121 && CANONDPDOFF(sourlo, 10) 122 && CANONDPDOFF(sourlo, 0)) return result; 123 #endif 124 } /* block */ 125 126 /* Loop to repair a non-canonical coefficent, as needed */ 127 inword=DECWORDS-1; /* current input word */ 128 uoff=0; /* bit offset of declet */ 129 encode=DFWORD(result, inword); 130 for (n=DECLETS-1; n>=0; n--) { /* count down declets of 10 bits */ 131 dpd=encode>>uoff; 132 uoff+=10; 133 if (uoff>32) { /* crossed uInt boundary */ 134 inword--; 135 encode=DFWORD(result, inword); 136 uoff-=32; 137 dpd|=encode<<(10-uoff); /* get pending bits */ 138 } 139 dpd&=0x3ff; /* clear uninteresting bits */ 140 if (dpd<0x16e) continue; /* must be canonical */ 141 canon=BIN2DPD[DPD2BIN[dpd]]; /* determine canonical declet */ 142 if (canon==dpd) continue; /* have canonical declet */ 143 /* need to replace declet */ 144 if (uoff>=10) { /* all within current word */ 145 encode&=~(0x3ff<<(uoff-10)); /* clear the 10 bits ready for replace */ 146 encode|=canon<<(uoff-10); /* insert the canonical form */ 147 DFWORD(result, inword)=encode; /* .. and save */ 148 continue; 149 } 150 /* straddled words */ 151 precode=DFWORD(result, inword+1); /* get previous */ 152 precode&=0xffffffff>>(10-uoff); /* clear top bits */ 153 DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff))); 154 encode&=0xffffffff<<uoff; /* clear bottom bits */ 155 encode|=canon>>(10-uoff); /* insert canonical */ 156 DFWORD(result, inword)=encode; /* .. and save */ 157 } /* n */ 158 return result; 159 } /* decCanonical */ 160 161/* ------------------------------------------------------------------ */ 162/* decDivide -- divide operations */ 163/* */ 164/* result gets the result of dividing dfl by dfr: */ 165/* dfl is the first decFloat (lhs) */ 166/* dfr is the second decFloat (rhs) */ 167/* set is the context */ 168/* op is the operation selector */ 169/* returns result */ 170/* */ 171/* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR. */ 172/* ------------------------------------------------------------------ */ 173#define DIVCOUNT 0 /* 1 to instrument subtractions counter */ 174#define DIVBASE ((uInt)BILLION) /* the base used for divide */ 175#define DIVOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */ 176#define DIVACCLEN (DIVOPLEN*3) /* accumulator length (ditto) */ 177static decFloat * decDivide(decFloat *result, const decFloat *dfl, 178 const decFloat *dfr, decContext *set, uInt op) { 179 decFloat quotient; /* for remainders */ 180 bcdnum num; /* for final conversion */ 181 uInt acc[DIVACCLEN]; /* coefficent in base-billion .. */ 182 uInt div[DIVOPLEN]; /* divisor in base-billion .. */ 183 uInt quo[DIVOPLEN+1]; /* quotient in base-billion .. */ 184 uByte bcdacc[(DIVOPLEN+1)*9+2]; /* for quotient in BCD, +1, +1 */ 185 uInt *msua, *msud, *msuq; /* -> msu of acc, div, and quo */ 186 Int divunits, accunits; /* lengths */ 187 Int quodigits; /* digits in quotient */ 188 uInt *lsua, *lsuq; /* -> current acc and quo lsus */ 189 Int length, multiplier; /* work */ 190 uInt carry, sign; /* .. */ 191 uInt *ua, *ud, *uq; /* .. */ 192 uByte *ub; /* .. */ 193 uInt uiwork; /* for macros */ 194 uInt divtop; /* top unit of div adjusted for estimating */ 195 #if DIVCOUNT 196 static uInt maxcount=0; /* worst-seen subtractions count */ 197 uInt divcount=0; /* subtractions count [this divide] */ 198 #endif 199 200 /* calculate sign */ 201 num.sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign; 202 203 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */ 204 /* NaNs are handled as usual */ 205 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set); 206 /* one or two infinities */ 207 if (DFISINF(dfl)) { 208 if (DFISINF(dfr)) return decInvalid(result, set); /* Two infinities bad */ 209 if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* as is rem */ 210 /* Infinity/x is infinite and quiet, even if x=0 */ 211 DFWORD(result, 0)=num.sign; 212 return decInfinity(result, result); 213 } 214 /* must be x/Infinity -- remainders are lhs */ 215 if (op&(REMAINDER|REMNEAR)) return decCanonical(result, dfl); 216 /* divides: return zero with correct sign and exponent depending */ 217 /* on op (Etiny for divide, 0 for divideInt) */ 218 decFloatZero(result); 219 if (op==DIVIDEINT) DFWORD(result, 0)|=num.sign; /* add sign */ 220 else DFWORD(result, 0)=num.sign; /* zeros the exponent, too */ 221 return result; 222 } 223 /* next, handle zero operands (x/0 and 0/x) */ 224 if (DFISZERO(dfr)) { /* x/0 */ 225 if (DFISZERO(dfl)) { /* 0/0 is undefined */ 226 decFloatZero(result); 227 DFWORD(result, 0)=DECFLOAT_qNaN; 228 set->status|=DEC_Division_undefined; 229 return result; 230 } 231 if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* bad rem */ 232 set->status|=DEC_Division_by_zero; 233 DFWORD(result, 0)=num.sign; 234 return decInfinity(result, result); /* x/0 -> signed Infinity */ 235 } 236 num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr); /* ideal exponent */ 237 if (DFISZERO(dfl)) { /* 0/x (x!=0) */ 238 /* if divide, result is 0 with ideal exponent; divideInt has */ 239 /* exponent=0, remainders give zero with lower exponent */ 240 if (op&DIVIDEINT) { 241 decFloatZero(result); 242 DFWORD(result, 0)|=num.sign; /* add sign */ 243 return result; 244 } 245 if (!(op&DIVIDE)) { /* a remainder */ 246 /* exponent is the minimum of the operands */ 247 num.exponent=MINI(GETEXPUN(dfl), GETEXPUN(dfr)); 248 /* if the result is zero the sign shall be sign of dfl */ 249 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign; 250 } 251 bcdacc[0]=0; 252 num.msd=bcdacc; /* -> 0 */ 253 num.lsd=bcdacc; /* .. */ 254 return decFinalize(result, &num, set); /* [divide may clamp exponent] */ 255 } /* 0/x */ 256 /* [here, both operands are known to be finite and non-zero] */ 257 258 /* extract the operand coefficents into 'units' which are */ 259 /* base-billion; the lhs is high-aligned in acc and the msu of both */ 260 /* acc and div is at the right-hand end of array (offset length-1); */ 261 /* the quotient can need one more unit than the operands as digits */ 262 /* in it are not necessarily aligned neatly; further, the quotient */ 263 /* may not start accumulating until after the end of the initial */ 264 /* operand in acc if that is small (e.g., 1) so the accumulator */ 265 /* must have at least that number of units extra (at the ls end) */ 266 GETCOEFFBILL(dfl, acc+DIVACCLEN-DIVOPLEN); 267 GETCOEFFBILL(dfr, div); 268 /* zero the low uInts of acc */ 269 acc[0]=0; 270 acc[1]=0; 271 acc[2]=0; 272 acc[3]=0; 273 #if DOUBLE 274 #if DIVOPLEN!=2 275 #error Unexpected Double DIVOPLEN 276 #endif 277 #elif QUAD 278 acc[4]=0; 279 acc[5]=0; 280 acc[6]=0; 281 acc[7]=0; 282 #if DIVOPLEN!=4 283 #error Unexpected Quad DIVOPLEN 284 #endif 285 #endif 286 287 /* set msu and lsu pointers */ 288 msua=acc+DIVACCLEN-1; /* [leading zeros removed below] */ 289 msuq=quo+DIVOPLEN; 290 /*[loop for div will terminate because operands are non-zero] */ 291 for (msud=div+DIVOPLEN-1; *msud==0;) msud--; 292 /* the initial least-significant unit of acc is set so acc appears */ 293 /* to have the same length as div. */ 294 /* This moves one position towards the least possible for each */ 295 /* iteration */ 296 divunits=(Int)(msud-div+1); /* precalculate */ 297 lsua=msua-divunits+1; /* initial working lsu of acc */ 298 lsuq=msuq; /* and of quo */ 299 300 /* set up the estimator for the multiplier; this is the msu of div, */ 301 /* plus two bits from the unit below (if any) rounded up by one if */ 302 /* there are any non-zero bits or units below that [the extra two */ 303 /* bits makes for a much better estimate when the top unit is small] */ 304 divtop=*msud<<2; 305 if (divunits>1) { 306 uInt *um=msud-1; 307 uInt d=*um; 308 if (d>=750000000) {divtop+=3; d-=750000000;} 309 else if (d>=500000000) {divtop+=2; d-=500000000;} 310 else if (d>=250000000) {divtop++; d-=250000000;} 311 if (d) divtop++; 312 else for (um--; um>=div; um--) if (*um) { 313 divtop++; 314 break; 315 } 316 } /* >1 unit */ 317 318 #if DECTRACE 319 {Int i; 320 printf("----- div="); 321 for (i=divunits-1; i>=0; i--) printf("%09ld ", (LI)div[i]); 322 printf("\n");} 323 #endif 324 325 /* now collect up to DECPMAX+1 digits in the quotient (this may */ 326 /* need OPLEN+1 uInts if unaligned) */ 327 quodigits=0; /* no digits yet */ 328 for (;; lsua--) { /* outer loop -- each input position */ 329 #if DECCHECK 330 if (lsua<acc) { 331 printf("Acc underrun...\n"); 332 break; 333 } 334 #endif 335 #if DECTRACE 336 printf("Outer: quodigits=%ld acc=", (LI)quodigits); 337 for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua); 338 printf("\n"); 339 #endif 340 *lsuq=0; /* default unit result is 0 */ 341 for (;;) { /* inner loop -- calculate quotient unit */ 342 /* strip leading zero units from acc (either there initially or */ 343 /* from subtraction below); this may strip all if exactly 0 */ 344 for (; *msua==0 && msua>=lsua;) msua--; 345 accunits=(Int)(msua-lsua+1); /* [maybe 0] */ 346 /* subtraction is only necessary and possible if there are as */ 347 /* least as many units remaining in acc for this iteration as */ 348 /* there are in div */ 349 if (accunits<divunits) { 350 if (accunits==0) msua++; /* restore */ 351 break; 352 } 353 354 /* If acc is longer than div then subtraction is definitely */ 355 /* possible (as msu of both is non-zero), but if they are the */ 356 /* same length a comparison is needed. */ 357 /* If a subtraction is needed then a good estimate of the */ 358 /* multiplier for the subtraction is also needed in order to */ 359 /* minimise the iterations of this inner loop because the */ 360 /* subtractions needed dominate division performance. */ 361 if (accunits==divunits) { 362 /* compare the high divunits of acc and div: */ 363 /* acc<div: this quotient unit is unchanged; subtraction */ 364 /* will be possible on the next iteration */ 365 /* acc==div: quotient gains 1, set acc=0 */ 366 /* acc>div: subtraction necessary at this position */ 367 for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break; 368 /* [now at first mismatch or lsu] */ 369 if (*ud>*ua) break; /* next time... */ 370 if (*ud==*ua) { /* all compared equal */ 371 *lsuq+=1; /* increment result */ 372 msua=lsua; /* collapse acc units */ 373 *msua=0; /* .. to a zero */ 374 break; 375 } 376 377 /* subtraction necessary; estimate multiplier [see above] */ 378 /* if both *msud and *msua are small it is cost-effective to */ 379 /* bring in part of the following units (if any) to get a */ 380 /* better estimate (assume some other non-zero in div) */ 381 #define DIVLO 1000000U 382 #define DIVHI (DIVBASE/DIVLO) 383 #if DECUSE64 384 if (divunits>1) { 385 /* there cannot be a *(msud-2) for DECDOUBLE so next is */ 386 /* an exact calculation unless DECQUAD (which needs to */ 387 /* assume bits out there if divunits>2) */ 388 uLong mul=(uLong)*msua * DIVBASE + *(msua-1); 389 uLong div=(uLong)*msud * DIVBASE + *(msud-1); 390 #if QUAD 391 if (divunits>2) div++; 392 #endif 393 mul/=div; 394 multiplier=(Int)mul; 395 } 396 else multiplier=*msua/(*msud); 397 #else 398 if (divunits>1 && *msua<DIVLO && *msud<DIVLO) { 399 multiplier=(*msua*DIVHI + *(msua-1)/DIVLO) 400 /(*msud*DIVHI + *(msud-1)/DIVLO +1); 401 } 402 else multiplier=(*msua<<2)/divtop; 403 #endif 404 } 405 else { /* accunits>divunits */ 406 /* msud is one unit 'lower' than msua, so estimate differently */ 407 #if DECUSE64 408 uLong mul; 409 /* as before, bring in extra digits if possible */ 410 if (divunits>1 && *msua<DIVLO && *msud<DIVLO) { 411 mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI 412 + *(msua-2)/DIVLO; 413 mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1); 414 } 415 else if (divunits==1) { 416 mul=(uLong)*msua * DIVBASE + *(msua-1); 417 mul/=*msud; /* no more to the right */ 418 } 419 else { 420 mul=(uLong)(*msua) * (uInt)(DIVBASE<<2) 421 + (*(msua-1)<<2); 422 mul/=divtop; /* [divtop already allows for sticky bits] */ 423 } 424 multiplier=(Int)mul; 425 #else 426 multiplier=*msua * ((DIVBASE<<2)/divtop); 427 #endif 428 } 429 if (multiplier==0) multiplier=1; /* marginal case */ 430 *lsuq+=multiplier; 431 432 #if DIVCOUNT 433 /* printf("Multiplier: %ld\n", (LI)multiplier); */ 434 divcount++; 435 #endif 436 437 /* Carry out the subtraction acc-(div*multiplier); for each */ 438 /* unit in div, do the multiply, split to units (see */ 439 /* decFloatMultiply for the algorithm), and subtract from acc */ 440 #define DIVMAGIC 2305843009U /* 2**61/10**9 */ 441 #define DIVSHIFTA 29 442 #define DIVSHIFTB 32 443 carry=0; 444 for (ud=div, ua=lsua; ud<=msud; ud++, ua++) { 445 uInt lo, hop; 446 #if DECUSE64 447 uLong sub=(uLong)multiplier*(*ud)+carry; 448 if (sub<DIVBASE) { 449 carry=0; 450 lo=(uInt)sub; 451 } 452 else { 453 hop=(uInt)(sub>>DIVSHIFTA); 454 carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB); 455 /* the estimate is now in hi; now calculate sub-hi*10**9 */ 456 /* to get the remainder (which will be <DIVBASE)) */ 457 lo=(uInt)sub; 458 lo-=carry*DIVBASE; /* low word of result */ 459 if (lo>=DIVBASE) { 460 lo-=DIVBASE; /* correct by +1 */ 461 carry++; 462 } 463 } 464 #else /* 32-bit */ 465 uInt hi; 466 /* calculate multiplier*(*ud) into hi and lo */ 467 LONGMUL32HI(hi, *ud, multiplier); /* get the high word */ 468 lo=multiplier*(*ud); /* .. and the low */ 469 lo+=carry; /* add the old hi */ 470 carry=hi+(lo<carry); /* .. with any carry */ 471 if (carry || lo>=DIVBASE) { /* split is needed */ 472 hop=(carry<<3)+(lo>>DIVSHIFTA); /* hi:lo/2**29 */ 473 LONGMUL32HI(carry, hop, DIVMAGIC); /* only need the high word */ 474 /* [DIVSHIFTB is 32, so carry can be used directly] */ 475 /* the estimate is now in carry; now calculate hi:lo-est*10**9; */ 476 /* happily the top word of the result is irrelevant because it */ 477 /* will always be zero so this needs only one multiplication */ 478 lo-=(carry*DIVBASE); 479 /* the correction here will be at most +1; do it */ 480 if (lo>=DIVBASE) { 481 lo-=DIVBASE; 482 carry++; 483 } 484 } 485 #endif 486 if (lo>*ua) { /* borrow needed */ 487 *ua+=DIVBASE; 488 carry++; 489 } 490 *ua-=lo; 491 } /* ud loop */ 492 if (carry) *ua-=carry; /* accdigits>divdigits [cannot borrow] */ 493 } /* inner loop */ 494 495 /* the outer loop terminates when there is either an exact result */ 496 /* or enough digits; first update the quotient digit count and */ 497 /* pointer (if any significant digits) */ 498 #if DECTRACE 499 if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq); 500 #endif 501 if (quodigits) { 502 quodigits+=9; /* had leading unit earlier */ 503 lsuq--; 504 if (quodigits>DECPMAX+1) break; /* have enough */ 505 } 506 else if (*lsuq) { /* first quotient digits */ 507 const uInt *pow; 508 for (pow=DECPOWERS; *lsuq>=*pow; pow++) quodigits++; 509 lsuq--; 510 /* [cannot have >DECPMAX+1 on first unit] */ 511 } 512 513 if (*msua!=0) continue; /* not an exact result */ 514 /* acc is zero iff used all of original units and zero down to lsua */ 515 /* (must also continue to original lsu for correct quotient length) */ 516 if (lsua>acc+DIVACCLEN-DIVOPLEN) continue; 517 for (; msua>lsua && *msua==0;) msua--; 518 if (*msua==0 && msua==lsua) break; 519 } /* outer loop */ 520 521 /* all of the original operand in acc has been covered at this point */ 522 /* quotient now has at least DECPMAX+2 digits */ 523 /* *msua is now non-0 if inexact and sticky bits */ 524 /* lsuq is one below the last uint of the quotient */ 525 lsuq++; /* set -> true lsu of quo */ 526 if (*msua) *lsuq|=1; /* apply sticky bit */ 527 528 /* quo now holds the (unrounded) quotient in base-billion; one */ 529 /* base-billion 'digit' per uInt. */ 530 #if DECTRACE 531 printf("DivQuo:"); 532 for (uq=msuq; uq>=lsuq; uq--) printf(" %09ld", (LI)*uq); 533 printf("\n"); 534 #endif 535 536 /* Now convert to BCD for rounding and cleanup, starting from the */ 537 /* most significant end [offset by one into bcdacc to leave room */ 538 /* for a possible carry digit if rounding for REMNEAR is needed] */ 539 for (uq=msuq, ub=bcdacc+1; uq>=lsuq; uq--, ub+=9) { 540 uInt top, mid, rem; /* work */ 541 if (*uq==0) { /* no split needed */ 542 UBFROMUI(ub, 0); /* clear 9 BCD8s */ 543 UBFROMUI(ub+4, 0); /* .. */ 544 *(ub+8)=0; /* .. */ 545 continue; 546 } 547 /* *uq is non-zero -- split the base-billion digit into */ 548 /* hi, mid, and low three-digits */ 549 #define divsplit9 1000000 /* divisor */ 550 #define divsplit6 1000 /* divisor */ 551 /* The splitting is done by simple divides and remainders, */ 552 /* assuming the compiler will optimize these [GCC does] */ 553 top=*uq/divsplit9; 554 rem=*uq%divsplit9; 555 mid=rem/divsplit6; 556 rem=rem%divsplit6; 557 /* lay out the nine BCD digits (plus one unwanted byte) */ 558 UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4])); 559 UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4])); 560 UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4])); 561 } /* BCD conversion loop */ 562 ub--; /* -> lsu */ 563 564 /* complete the bcdnum; quodigits is correct, so the position of */ 565 /* the first non-zero is known */ 566 num.msd=bcdacc+1+(msuq-lsuq+1)*9-quodigits; 567 num.lsd=ub; 568 569 /* make exponent adjustments, etc */ 570 if (lsua<acc+DIVACCLEN-DIVOPLEN) { /* used extra digits */ 571 num.exponent-=(Int)((acc+DIVACCLEN-DIVOPLEN-lsua)*9); 572 /* if the result was exact then there may be up to 8 extra */ 573 /* trailing zeros in the overflowed quotient final unit */ 574 if (*msua==0) { 575 for (; *ub==0;) ub--; /* drop zeros */ 576 num.exponent+=(Int)(num.lsd-ub); /* and adjust exponent */ 577 num.lsd=ub; 578 } 579 } /* adjustment needed */ 580 581 #if DIVCOUNT 582 if (divcount>maxcount) { /* new high-water nark */ 583 maxcount=divcount; 584 printf("DivNewMaxCount: %ld\n", (LI)maxcount); 585 } 586 #endif 587 588 if (op&DIVIDE) return decFinalize(result, &num, set); /* all done */ 589 590 /* Is DIVIDEINT or a remainder; there is more to do -- first form */ 591 /* the integer (this is done 'after the fact', unlike as in */ 592 /* decNumber, so as not to tax DIVIDE) */ 593 594 /* The first non-zero digit will be in the first 9 digits, known */ 595 /* from quodigits and num.msd, so there is always space for DECPMAX */ 596 /* digits */ 597 598 length=(Int)(num.lsd-num.msd+1); 599 /*printf("Length exp: %ld %ld\n", (LI)length, (LI)num.exponent); */ 600 601 if (length+num.exponent>DECPMAX) { /* cannot fit */ 602 decFloatZero(result); 603 DFWORD(result, 0)=DECFLOAT_qNaN; 604 set->status|=DEC_Division_impossible; 605 return result; 606 } 607 608 if (num.exponent>=0) { /* already an int, or need pad zeros */ 609 for (ub=num.lsd+1; ub<=num.lsd+num.exponent; ub++) *ub=0; 610 num.lsd+=num.exponent; 611 } 612 else { /* too long: round or truncate needed */ 613 Int drop=-num.exponent; 614 if (!(op&REMNEAR)) { /* simple truncate */ 615 num.lsd-=drop; 616 if (num.lsd<num.msd) { /* truncated all */ 617 num.lsd=num.msd; /* make 0 */ 618 *num.lsd=0; /* .. [sign still relevant] */ 619 } 620 } 621 else { /* round to nearest even [sigh] */ 622 /* round-to-nearest, in-place; msd is at or to right of bcdacc+1 */ 623 /* (this is a special case of Quantize -- q.v. for commentary) */ 624 uByte *roundat; /* -> re-round digit */ 625 uByte reround; /* reround value */ 626 *(num.msd-1)=0; /* in case of left carry, or make 0 */ 627 if (drop<length) roundat=num.lsd-drop+1; 628 else if (drop==length) roundat=num.msd; 629 else roundat=num.msd-1; /* [-> 0] */ 630 reround=*roundat; 631 for (ub=roundat+1; ub<=num.lsd; ub++) { 632 if (*ub!=0) { 633 reround=DECSTICKYTAB[reround]; 634 break; 635 } 636 } /* check stickies */ 637 if (roundat>num.msd) num.lsd=roundat-1; 638 else { 639 num.msd--; /* use the 0 .. */ 640 num.lsd=num.msd; /* .. at the new MSD place */ 641 } 642 if (reround!=0) { /* discarding non-zero */ 643 uInt bump=0; 644 /* rounding is DEC_ROUND_HALF_EVEN always */ 645 if (reround>5) bump=1; /* >0.5 goes up */ 646 else if (reround==5) /* exactly 0.5000 .. */ 647 bump=*(num.lsd) & 0x01; /* .. up iff [new] lsd is odd */ 648 if (bump!=0) { /* need increment */ 649 /* increment the coefficient; this might end up with 1000... */ 650 ub=num.lsd; 651 for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0); 652 for (; *ub==9; ub--) *ub=0; /* at most 3 more */ 653 *ub+=1; 654 if (ub<num.msd) num.msd--; /* carried */ 655 } /* bump needed */ 656 } /* reround!=0 */ 657 } /* remnear */ 658 } /* round or truncate needed */ 659 num.exponent=0; /* all paths */ 660 /*decShowNum(&num, "int"); */ 661 662 if (op&DIVIDEINT) return decFinalize(result, &num, set); /* all done */ 663 664 /* Have a remainder to calculate */ 665 decFinalize("ient, &num, set); /* lay out the integer so far */ 666 DFWORD("ient, 0)^=DECFLOAT_Sign; /* negate it */ 667 sign=DFWORD(dfl, 0); /* save sign of dfl */ 668 decFloatFMA(result, "ient, dfr, dfl, set); 669 if (!DFISZERO(result)) return result; 670 /* if the result is zero the sign shall be sign of dfl */ 671 DFWORD("ient, 0)=sign; /* construct decFloat of sign */ 672 return decFloatCopySign(result, result, "ient); 673 } /* decDivide */ 674 675/* ------------------------------------------------------------------ */ 676/* decFiniteMultiply -- multiply two finite decFloats */ 677/* */ 678/* num gets the result of multiplying dfl and dfr */ 679/* bcdacc .. with the coefficient in this array */ 680/* dfl is the first decFloat (lhs) */ 681/* dfr is the second decFloat (rhs) */ 682/* */ 683/* This effects the multiplication of two decFloats, both known to be */ 684/* finite, leaving the result in a bcdnum ready for decFinalize (for */ 685/* use in Multiply) or in a following addition (FMA). */ 686/* */ 687/* bcdacc must have space for at least DECPMAX9*18+1 bytes. */ 688/* No error is possible and no status is set. */ 689/* ------------------------------------------------------------------ */ 690/* This routine has two separate implementations of the core */ 691/* multiplication; both using base-billion. One uses only 32-bit */ 692/* variables (Ints and uInts) or smaller; the other uses uLongs (for */ 693/* multiplication and addition only). Both implementations cover */ 694/* both arithmetic sizes (DOUBLE and QUAD) in order to allow timing */ 695/* comparisons. In any one compilation only one implementation for */ 696/* each size can be used, and if DECUSE64 is 0 then use of the 32-bit */ 697/* version is forced. */ 698/* */ 699/* Historical note: an earlier version of this code also supported the */ 700/* 256-bit format and has been preserved. That is somewhat trickier */ 701/* during lazy carry splitting because the initial quotient estimate */ 702/* (est) can exceed 32 bits. */ 703 704#define MULTBASE ((uInt)BILLION) /* the base used for multiply */ 705#define MULOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */ 706#define MULACCLEN (MULOPLEN*2) /* accumulator length (ditto) */ 707#define LEADZEROS (MULACCLEN*9 - DECPMAX*2) /* leading zeros always */ 708 709/* Assertions: exponent not too large and MULACCLEN is a multiple of 4 */ 710#if DECEMAXD>9 711 #error Exponent may overflow when doubled for Multiply 712#endif 713#if MULACCLEN!=(MULACCLEN/4)*4 714 /* This assumption is used below only for initialization */ 715 #error MULACCLEN is not a multiple of 4 716#endif 717 718static void decFiniteMultiply(bcdnum *num, uByte *bcdacc, 719 const decFloat *dfl, const decFloat *dfr) { 720 uInt bufl[MULOPLEN]; /* left coefficient (base-billion) */ 721 uInt bufr[MULOPLEN]; /* right coefficient (base-billion) */ 722 uInt *ui, *uj; /* work */ 723 uByte *ub; /* .. */ 724 uInt uiwork; /* for macros */ 725 726 #if DECUSE64 727 uLong accl[MULACCLEN]; /* lazy accumulator (base-billion+) */ 728 uLong *pl; /* work -> lazy accumulator */ 729 uInt acc[MULACCLEN]; /* coefficent in base-billion .. */ 730 #else 731 uInt acc[MULACCLEN*2]; /* accumulator in base-billion .. */ 732 #endif 733 uInt *pa; /* work -> accumulator */ 734 /*printf("Base10**9: OpLen=%d MulAcclen=%d\n", OPLEN, MULACCLEN); */ 735 736 /* Calculate sign and exponent */ 737 num->sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign; 738 num->exponent=GETEXPUN(dfl)+GETEXPUN(dfr); /* [see assertion above] */ 739 740 /* Extract the coefficients and prepare the accumulator */ 741 /* the coefficients of the operands are decoded into base-billion */ 742 /* numbers in uInt arrays (bufl and bufr, LSD at offset 0) of the */ 743 /* appropriate size. */ 744 GETCOEFFBILL(dfl, bufl); 745 GETCOEFFBILL(dfr, bufr); 746 #if DECTRACE && 0 747 printf("CoeffbL:"); 748 for (ui=bufl+MULOPLEN-1; ui>=bufl; ui--) printf(" %08lx", (LI)*ui); 749 printf("\n"); 750 printf("CoeffbR:"); 751 for (uj=bufr+MULOPLEN-1; uj>=bufr; uj--) printf(" %08lx", (LI)*uj); 752 printf("\n"); 753 #endif 754 755 /* start the 64-bit/32-bit differing paths... */ 756#if DECUSE64 757 758 /* zero the accumulator */ 759 #if MULACCLEN==4 760 accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0; 761 #else /* use a loop */ 762 /* MULACCLEN is a multiple of four, asserted above */ 763 for (pl=accl; pl<accl+MULACCLEN; pl+=4) { 764 *pl=0; *(pl+1)=0; *(pl+2)=0; *(pl+3)=0;/* [reduce overhead] */ 765 } /* pl */ 766 #endif 767 768 /* Effect the multiplication */ 769 /* The multiplcation proceeds using MFC's lazy-carry resolution */ 770 /* algorithm from decNumber. First, the multiplication is */ 771 /* effected, allowing accumulation of the partial products (which */ 772 /* are in base-billion at each column position) into 64 bits */ 773 /* without resolving back to base=billion after each addition. */ 774 /* These 64-bit numbers (which may contain up to 19 decimal digits) */ 775 /* are then split using the Clark & Cowlishaw algorithm (see below). */ 776 /* [Testing for 0 in the inner loop is not really a 'win'] */ 777 for (ui=bufr; ui<bufr+MULOPLEN; ui++) { /* over each item in rhs */ 778 if (*ui==0) continue; /* product cannot affect result */ 779 pl=accl+(ui-bufr); /* where to add the lhs */ 780 for (uj=bufl; uj<bufl+MULOPLEN; uj++, pl++) { /* over each item in lhs */ 781 /* if (*uj==0) continue; // product cannot affect result */ 782 *pl+=((uLong)*ui)*(*uj); 783 } /* uj */ 784 } /* ui */ 785 786 /* The 64-bit carries must now be resolved; this means that a */ 787 /* quotient/remainder has to be calculated for base-billion (1E+9). */ 788 /* For this, Clark & Cowlishaw's quotient estimation approach (also */ 789 /* used in decNumber) is needed, because 64-bit divide is generally */ 790 /* extremely slow on 32-bit machines, and may be slower than this */ 791 /* approach even on 64-bit machines. This algorithm splits X */ 792 /* using: */ 793 /* */ 794 /* magic=2**(A+B)/1E+9; // 'magic number' */ 795 /* hop=X/2**A; // high order part of X (by shift) */ 796 /* est=magic*hop/2**B // quotient estimate (may be low by 1) */ 797 /* */ 798 /* A and B are quite constrained; hop and magic must fit in 32 bits, */ 799 /* and 2**(A+B) must be as large as possible (which is 2**61 if */ 800 /* magic is to fit). Further, maxX increases with the length of */ 801 /* the operands (and hence the number of partial products */ 802 /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */ 803 /* */ 804 /* It can be shown that when OPLEN is 2 then the maximum error in */ 805 /* the estimated quotient is <1, but for larger maximum x the */ 806 /* maximum error is above 1 so a correction that is >1 may be */ 807 /* needed. Values of A and B are chosen to satisfy the constraints */ 808 /* just mentioned while minimizing the maximum error (and hence the */ 809 /* maximum correction), as shown in the following table: */ 810 /* */ 811 /* Type OPLEN A B maxX maxError maxCorrection */ 812 /* --------------------------------------------------------- */ 813 /* DOUBLE 2 29 32 <2*10**18 0.63 1 */ 814 /* QUAD 4 30 31 <4*10**18 1.17 2 */ 815 /* */ 816 /* In the OPLEN==2 case there is most choice, but the value for B */ 817 /* of 32 has a big advantage as then the calculation of the */ 818 /* estimate requires no shifting; the compiler can extract the high */ 819 /* word directly after multiplying magic*hop. */ 820 #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */ 821 #if DOUBLE 822 #define MULSHIFTA 29 823 #define MULSHIFTB 32 824 #elif QUAD 825 #define MULSHIFTA 30 826 #define MULSHIFTB 31 827 #else 828 #error Unexpected type 829 #endif 830 831 #if DECTRACE 832 printf("MulAccl:"); 833 for (pl=accl+MULACCLEN-1; pl>=accl; pl--) 834 printf(" %08lx:%08lx", (LI)(*pl>>32), (LI)(*pl&0xffffffff)); 835 printf("\n"); 836 #endif 837 838 for (pl=accl, pa=acc; pl<accl+MULACCLEN; pl++, pa++) { /* each column position */ 839 uInt lo, hop; /* work */ 840 uInt est; /* cannot exceed 4E+9 */ 841 if (*pl>=MULTBASE) { 842 /* *pl holds a binary number which needs to be split */ 843 hop=(uInt)(*pl>>MULSHIFTA); 844 est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB); 845 /* the estimate is now in est; now calculate hi:lo-est*10**9; */ 846 /* happily the top word of the result is irrelevant because it */ 847 /* will always be zero so this needs only one multiplication */ 848 lo=(uInt)(*pl-((uLong)est*MULTBASE)); /* low word of result */ 849 /* If QUAD, the correction here could be +2 */ 850 if (lo>=MULTBASE) { 851 lo-=MULTBASE; /* correct by +1 */ 852 est++; 853 #if QUAD 854 /* may need to correct by +2 */ 855 if (lo>=MULTBASE) { 856 lo-=MULTBASE; 857 est++; 858 } 859 #endif 860 } 861 /* finally place lo as the new coefficient 'digit' and add est to */ 862 /* the next place up [this is safe because this path is never */ 863 /* taken on the final iteration as *pl will fit] */ 864 *pa=lo; 865 *(pl+1)+=est; 866 } /* *pl needed split */ 867 else { /* *pl<MULTBASE */ 868 *pa=(uInt)*pl; /* just copy across */ 869 } 870 } /* pl loop */ 871 872#else /* 32-bit */ 873 for (pa=acc;; pa+=4) { /* zero the accumulator */ 874 *pa=0; *(pa+1)=0; *(pa+2)=0; *(pa+3)=0; /* [reduce overhead] */ 875 if (pa==acc+MULACCLEN*2-4) break; /* multiple of 4 asserted */ 876 } /* pa */ 877 878 /* Effect the multiplication */ 879 /* uLongs are not available (and in particular, there is no uLong */ 880 /* divide) but it is still possible to use MFC's lazy-carry */ 881 /* resolution algorithm from decNumber. First, the multiplication */ 882 /* is effected, allowing accumulation of the partial products */ 883 /* (which are in base-billion at each column position) into 64 bits */ 884 /* [with the high-order 32 bits in each position being held at */ 885 /* offset +ACCLEN from the low-order 32 bits in the accumulator]. */ 886 /* These 64-bit numbers (which may contain up to 19 decimal digits) */ 887 /* are then split using the Clark & Cowlishaw algorithm (see */ 888 /* below). */ 889 for (ui=bufr;; ui++) { /* over each item in rhs */ 890 uInt hi, lo; /* words of exact multiply result */ 891 pa=acc+(ui-bufr); /* where to add the lhs */ 892 for (uj=bufl;; uj++, pa++) { /* over each item in lhs */ 893 LONGMUL32HI(hi, *ui, *uj); /* calculate product of digits */ 894 lo=(*ui)*(*uj); /* .. */ 895 *pa+=lo; /* accumulate low bits and .. */ 896 *(pa+MULACCLEN)+=hi+(*pa<lo); /* .. high bits with any carry */ 897 if (uj==bufl+MULOPLEN-1) break; 898 } 899 if (ui==bufr+MULOPLEN-1) break; 900 } 901 902 /* The 64-bit carries must now be resolved; this means that a */ 903 /* quotient/remainder has to be calculated for base-billion (1E+9). */ 904 /* For this, Clark & Cowlishaw's quotient estimation approach (also */ 905 /* used in decNumber) is needed, because 64-bit divide is generally */ 906 /* extremely slow on 32-bit machines. This algorithm splits X */ 907 /* using: */ 908 /* */ 909 /* magic=2**(A+B)/1E+9; // 'magic number' */ 910 /* hop=X/2**A; // high order part of X (by shift) */ 911 /* est=magic*hop/2**B // quotient estimate (may be low by 1) */ 912 /* */ 913 /* A and B are quite constrained; hop and magic must fit in 32 bits, */ 914 /* and 2**(A+B) must be as large as possible (which is 2**61 if */ 915 /* magic is to fit). Further, maxX increases with the length of */ 916 /* the operands (and hence the number of partial products */ 917 /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */ 918 /* */ 919 /* It can be shown that when OPLEN is 2 then the maximum error in */ 920 /* the estimated quotient is <1, but for larger maximum x the */ 921 /* maximum error is above 1 so a correction that is >1 may be */ 922 /* needed. Values of A and B are chosen to satisfy the constraints */ 923 /* just mentioned while minimizing the maximum error (and hence the */ 924 /* maximum correction), as shown in the following table: */ 925 /* */ 926 /* Type OPLEN A B maxX maxError maxCorrection */ 927 /* --------------------------------------------------------- */ 928 /* DOUBLE 2 29 32 <2*10**18 0.63 1 */ 929 /* QUAD 4 30 31 <4*10**18 1.17 2 */ 930 /* */ 931 /* In the OPLEN==2 case there is most choice, but the value for B */ 932 /* of 32 has a big advantage as then the calculation of the */ 933 /* estimate requires no shifting; the high word is simply */ 934 /* calculated from multiplying magic*hop. */ 935 #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */ 936 #if DOUBLE 937 #define MULSHIFTA 29 938 #define MULSHIFTB 32 939 #elif QUAD 940 #define MULSHIFTA 30 941 #define MULSHIFTB 31 942 #else 943 #error Unexpected type 944 #endif 945 946 #if DECTRACE 947 printf("MulHiLo:"); 948 for (pa=acc+MULACCLEN-1; pa>=acc; pa--) 949 printf(" %08lx:%08lx", (LI)*(pa+MULACCLEN), (LI)*pa); 950 printf("\n"); 951 #endif 952 953 for (pa=acc;; pa++) { /* each low uInt */ 954 uInt hi, lo; /* words of exact multiply result */ 955 uInt hop, estlo; /* work */ 956 #if QUAD 957 uInt esthi; /* .. */ 958 #endif 959 960 lo=*pa; 961 hi=*(pa+MULACCLEN); /* top 32 bits */ 962 /* hi and lo now hold a binary number which needs to be split */ 963 964 #if DOUBLE 965 hop=(hi<<3)+(lo>>MULSHIFTA); /* hi:lo/2**29 */ 966 LONGMUL32HI(estlo, hop, MULMAGIC);/* only need the high word */ 967 /* [MULSHIFTB is 32, so estlo can be used directly] */ 968 /* the estimate is now in estlo; now calculate hi:lo-est*10**9; */ 969 /* happily the top word of the result is irrelevant because it */ 970 /* will always be zero so this needs only one multiplication */ 971 lo-=(estlo*MULTBASE); 972 /* esthi=0; // high word is ignored below */ 973 /* the correction here will be at most +1; do it */ 974 if (lo>=MULTBASE) { 975 lo-=MULTBASE; 976 estlo++; 977 } 978 #elif QUAD 979 hop=(hi<<2)+(lo>>MULSHIFTA); /* hi:lo/2**30 */ 980 LONGMUL32HI(esthi, hop, MULMAGIC);/* shift will be 31 .. */ 981 estlo=hop*MULMAGIC; /* .. so low word needed */ 982 estlo=(esthi<<1)+(estlo>>MULSHIFTB); /* [just the top bit] */ 983 /* esthi=0; // high word is ignored below */ 984 lo-=(estlo*MULTBASE); /* as above */ 985 /* the correction here could be +1 or +2 */ 986 if (lo>=MULTBASE) { 987 lo-=MULTBASE; 988 estlo++; 989 } 990 if (lo>=MULTBASE) { 991 lo-=MULTBASE; 992 estlo++; 993 } 994 #else 995 #error Unexpected type 996 #endif 997 998 /* finally place lo as the new accumulator digit and add est to */ 999 /* the next place up; this latter add could cause a carry of 1 */ 1000 /* to the high word of the next place */ 1001 *pa=lo; 1002 *(pa+1)+=estlo; 1003 /* esthi is always 0 for DOUBLE and QUAD so this is skipped */ 1004 /* *(pa+1+MULACCLEN)+=esthi; */ 1005 if (*(pa+1)<estlo) *(pa+1+MULACCLEN)+=1; /* carry */ 1006 if (pa==acc+MULACCLEN-2) break; /* [MULACCLEN-1 will never need split] */ 1007 } /* pa loop */ 1008#endif 1009 1010 /* At this point, whether using the 64-bit or the 32-bit paths, the */ 1011 /* accumulator now holds the (unrounded) result in base-billion; */ 1012 /* one base-billion 'digit' per uInt. */ 1013 #if DECTRACE 1014 printf("MultAcc:"); 1015 for (pa=acc+MULACCLEN-1; pa>=acc; pa--) printf(" %09ld", (LI)*pa); 1016 printf("\n"); 1017 #endif 1018 1019 /* Now convert to BCD for rounding and cleanup, starting from the */ 1020 /* most significant end */ 1021 pa=acc+MULACCLEN-1; 1022 if (*pa!=0) num->msd=bcdacc+LEADZEROS;/* drop known lead zeros */ 1023 else { /* >=1 word of leading zeros */ 1024 num->msd=bcdacc; /* known leading zeros are gone */ 1025 pa--; /* skip first word .. */ 1026 for (; *pa==0; pa--) if (pa==acc) break; /* .. and any more leading 0s */ 1027 } 1028 for (ub=bcdacc;; pa--, ub+=9) { 1029 if (*pa!=0) { /* split(s) needed */ 1030 uInt top, mid, rem; /* work */ 1031 /* *pa is non-zero -- split the base-billion acc digit into */ 1032 /* hi, mid, and low three-digits */ 1033 #define mulsplit9 1000000 /* divisor */ 1034 #define mulsplit6 1000 /* divisor */ 1035 /* The splitting is done by simple divides and remainders, */ 1036 /* assuming the compiler will optimize these where useful */ 1037 /* [GCC does] */ 1038 top=*pa/mulsplit9; 1039 rem=*pa%mulsplit9; 1040 mid=rem/mulsplit6; 1041 rem=rem%mulsplit6; 1042 /* lay out the nine BCD digits (plus one unwanted byte) */ 1043 UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4])); 1044 UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4])); 1045 UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4])); 1046 } 1047 else { /* *pa==0 */ 1048 UBFROMUI(ub, 0); /* clear 9 BCD8s */ 1049 UBFROMUI(ub+4, 0); /* .. */ 1050 *(ub+8)=0; /* .. */ 1051 } 1052 if (pa==acc) break; 1053 } /* BCD conversion loop */ 1054 1055 num->lsd=ub+8; /* complete the bcdnum .. */ 1056 1057 #if DECTRACE 1058 decShowNum(num, "postmult"); 1059 decFloatShow(dfl, "dfl"); 1060 decFloatShow(dfr, "dfr"); 1061 #endif 1062 return; 1063 } /* decFiniteMultiply */ 1064 1065/* ------------------------------------------------------------------ */ 1066/* decFloatAbs -- absolute value, heeding NaNs, etc. */ 1067/* */ 1068/* result gets the canonicalized df with sign 0 */ 1069/* df is the decFloat to abs */ 1070/* set is the context */ 1071/* returns result */ 1072/* */ 1073/* This has the same effect as decFloatPlus unless df is negative, */ 1074/* in which case it has the same effect as decFloatMinus. The */ 1075/* effect is also the same as decFloatCopyAbs except that NaNs are */ 1076/* handled normally (the sign of a NaN is not affected, and an sNaN */ 1077/* will signal) and the result will be canonical. */ 1078/* ------------------------------------------------------------------ */ 1079decFloat * decFloatAbs(decFloat *result, const decFloat *df, 1080 decContext *set) { 1081 if (DFISNAN(df)) return decNaNs(result, df, NULL, set); 1082 decCanonical(result, df); /* copy and check */ 1083 DFBYTE(result, 0)&=~0x80; /* zero sign bit */ 1084 return result; 1085 } /* decFloatAbs */ 1086 1087/* ------------------------------------------------------------------ */ 1088/* decFloatAdd -- add two decFloats */ 1089/* */ 1090/* result gets the result of adding dfl and dfr: */ 1091/* dfl is the first decFloat (lhs) */ 1092/* dfr is the second decFloat (rhs) */ 1093/* set is the context */ 1094/* returns result */ 1095/* */ 1096/* ------------------------------------------------------------------ */ 1097#if QUAD 1098/* Table for testing MSDs for fastpath elimination; returns the MSD of */ 1099/* a decDouble or decQuad (top 6 bits tested) ignoring the sign. */ 1100/* Infinities return -32 and NaNs return -128 so that summing the two */ 1101/* MSDs also allows rapid tests for the Specials (see code below). */ 1102const Int DECTESTMSD[64]={ 1103 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 1104 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128, 1105 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 1106 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128}; 1107#else 1108/* The table for testing MSDs is shared between the modules */ 1109extern const Int DECTESTMSD[64]; 1110#endif 1111 1112decFloat * decFloatAdd(decFloat *result, 1113 const decFloat *dfl, const decFloat *dfr, 1114 decContext *set) { 1115 bcdnum num; /* for final conversion */ 1116 Int bexpl, bexpr; /* left and right biased exponents */ 1117 uByte *ub, *us, *ut; /* work */ 1118 uInt uiwork; /* for macros */ 1119 #if QUAD 1120 uShort uswork; /* .. */ 1121 #endif 1122 1123 uInt sourhil, sourhir; /* top words from source decFloats */ 1124 /* [valid only through end of */ 1125 /* fastpath code -- before swap] */ 1126 uInt diffsign; /* non-zero if signs differ */ 1127 uInt carry; /* carry: 0 or 1 before add loop */ 1128 Int overlap; /* coefficient overlap (if full) */ 1129 Int summ; /* sum of the MSDs */ 1130 /* the following buffers hold coefficients with various alignments */ 1131 /* (see commentary and diagrams below) */ 1132 uByte acc[4+2+DECPMAX*3+8]; 1133 uByte buf[4+2+DECPMAX*2]; 1134 uByte *umsd, *ulsd; /* local MSD and LSD pointers */ 1135 1136 #if DECLITEND 1137 #define CARRYPAT 0x01000000 /* carry=1 pattern */ 1138 #else 1139 #define CARRYPAT 0x00000001 /* carry=1 pattern */ 1140 #endif 1141 1142 /* Start decoding the arguments */ 1143 /* The initial exponents are placed into the opposite Ints to */ 1144 /* that which might be expected; there are two sets of data to */ 1145 /* keep track of (each decFloat and the corresponding exponent), */ 1146 /* and this scheme means that at the swap point (after comparing */ 1147 /* exponents) only one pair of words needs to be swapped */ 1148 /* whichever path is taken (thereby minimising worst-case path). */ 1149 /* The calculated exponents will be nonsense when the arguments are */ 1150 /* Special, but are not used in that path */ 1151 sourhil=DFWORD(dfl, 0); /* LHS top word */ 1152 summ=DECTESTMSD[sourhil>>26]; /* get first MSD for testing */ 1153 bexpr=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */ 1154 bexpr+=GETECON(dfl); /* .. + continuation */ 1155 1156 sourhir=DFWORD(dfr, 0); /* RHS top word */ 1157 summ+=DECTESTMSD[sourhir>>26]; /* sum MSDs for testing */ 1158 bexpl=DECCOMBEXP[sourhir>>26]; 1159 bexpl+=GETECON(dfr); 1160 1161 /* here bexpr has biased exponent from lhs, and vice versa */ 1162 1163 diffsign=(sourhil^sourhir)&DECFLOAT_Sign; 1164 1165 /* now determine whether to take a fast path or the full-function */ 1166 /* slow path. The slow path must be taken when: */ 1167 /* -- both numbers are finite, and: */ 1168 /* the exponents are different, or */ 1169 /* the signs are different, or */ 1170 /* the sum of the MSDs is >8 (hence might overflow) */ 1171 /* specialness and the sum of the MSDs can be tested at once using */ 1172 /* the summ value just calculated, so the test for specials is no */ 1173 /* longer on the worst-case path (as of 3.60) */ 1174 1175 if (summ<=8) { /* MSD+MSD is good, or there is a special */ 1176 if (summ<0) { /* there is a special */ 1177 /* Inf+Inf would give -64; Inf+finite is -32 or higher */ 1178 if (summ<-64) return decNaNs(result, dfl, dfr, set); /* one or two NaNs */ 1179 /* two infinities with different signs is invalid */ 1180 if (summ==-64 && diffsign) return decInvalid(result, set); 1181 if (DFISINF(dfl)) return decInfinity(result, dfl); /* LHS is infinite */ 1182 return decInfinity(result, dfr); /* RHS must be Inf */ 1183 } 1184 /* Here when both arguments are finite; fast path is possible */ 1185 /* (currently only for aligned and same-sign) */ 1186 if (bexpr==bexpl && !diffsign) { 1187 uInt tac[DECLETS+1]; /* base-1000 coefficient */ 1188 uInt encode; /* work */ 1189 1190 /* Get one coefficient as base-1000 and add the other */ 1191 GETCOEFFTHOU(dfl, tac); /* least-significant goes to [0] */ 1192 ADDCOEFFTHOU(dfr, tac); 1193 /* here the sum of the MSDs (plus any carry) will be <10 due to */ 1194 /* the fastpath test earlier */ 1195 1196 /* construct the result; low word is the same for both formats */ 1197 encode =BIN2DPD[tac[0]]; 1198 encode|=BIN2DPD[tac[1]]<<10; 1199 encode|=BIN2DPD[tac[2]]<<20; 1200 encode|=BIN2DPD[tac[3]]<<30; 1201 DFWORD(result, (DECBYTES/4)-1)=encode; 1202 1203 /* collect next two declets (all that remains, for Double) */ 1204 encode =BIN2DPD[tac[3]]>>2; 1205 encode|=BIN2DPD[tac[4]]<<8; 1206 1207 #if QUAD 1208 /* complete and lay out middling words */ 1209 encode|=BIN2DPD[tac[5]]<<18; 1210 encode|=BIN2DPD[tac[6]]<<28; 1211 DFWORD(result, 2)=encode; 1212 1213 encode =BIN2DPD[tac[6]]>>4; 1214 encode|=BIN2DPD[tac[7]]<<6; 1215 encode|=BIN2DPD[tac[8]]<<16; 1216 encode|=BIN2DPD[tac[9]]<<26; 1217 DFWORD(result, 1)=encode; 1218 1219 /* and final two declets */ 1220 encode =BIN2DPD[tac[9]]>>6; 1221 encode|=BIN2DPD[tac[10]]<<4; 1222 #endif 1223 1224 /* add exponent continuation and sign (from either argument) */ 1225 encode|=sourhil & (ECONMASK | DECFLOAT_Sign); 1226 1227 /* create lookup index = MSD + top two bits of biased exponent <<4 */ 1228 tac[DECLETS]|=(bexpl>>DECECONL)<<4; 1229 encode|=DECCOMBFROM[tac[DECLETS]]; /* add constructed combination field */ 1230 DFWORD(result, 0)=encode; /* complete */ 1231 1232 /* decFloatShow(result, ">"); */ 1233 return result; 1234 } /* fast path OK */ 1235 /* drop through to slow path */ 1236 } /* low sum or Special(s) */ 1237 1238 /* Slow path required -- arguments are finite and might overflow, */ 1239 /* or require alignment, or might have different signs */ 1240 1241 /* now swap either exponents or argument pointers */ 1242 if (bexpl<=bexpr) { 1243 /* original left is bigger */ 1244 Int bexpswap=bexpl; 1245 bexpl=bexpr; 1246 bexpr=bexpswap; 1247 /* printf("left bigger\n"); */ 1248 } 1249 else { 1250 const decFloat *dfswap=dfl; 1251 dfl=dfr; 1252 dfr=dfswap; 1253 /* printf("right bigger\n"); */ 1254 } 1255 /* [here dfl and bexpl refer to the datum with the larger exponent, */ 1256 /* of if the exponents are equal then the original LHS argument] */ 1257 1258 /* if lhs is zero then result will be the rhs (now known to have */ 1259 /* the smaller exponent), which also may need to be tested for zero */ 1260 /* for the weird IEEE 754 sign rules */ 1261 if (DFISZERO(dfl)) { 1262 decCanonical(result, dfr); /* clean copy */ 1263 /* "When the sum of two operands with opposite signs is */ 1264 /* exactly zero, the sign of that sum shall be '+' in all */ 1265 /* rounding modes except round toward -Infinity, in which */ 1266 /* mode that sign shall be '-'." */ 1267 if (diffsign && DFISZERO(result)) { 1268 DFWORD(result, 0)&=~DECFLOAT_Sign; /* assume sign 0 */ 1269 if (set->round==DEC_ROUND_FLOOR) DFWORD(result, 0)|=DECFLOAT_Sign; 1270 } 1271 return result; 1272 } /* numfl is zero */ 1273 /* [here, LHS is non-zero; code below assumes that] */ 1274 1275 /* Coefficients layout during the calculations to follow: */ 1276 /* */ 1277 /* Overlap case: */ 1278 /* +------------------------------------------------+ */ 1279 /* acc: |0000| coeffa | tail B | | */ 1280 /* +------------------------------------------------+ */ 1281 /* buf: |0000| pad0s | coeffb | | */ 1282 /* +------------------------------------------------+ */ 1283 /* */ 1284 /* Touching coefficients or gap: */ 1285 /* +------------------------------------------------+ */ 1286 /* acc: |0000| coeffa | gap | coeffb | */ 1287 /* +------------------------------------------------+ */ 1288 /* [buf not used or needed; gap clamped to Pmax] */ 1289 1290 /* lay out lhs coefficient into accumulator; this starts at acc+4 */ 1291 /* for decDouble or acc+6 for decQuad so the LSD is word- */ 1292 /* aligned; the top word gap is there only in case a carry digit */ 1293 /* is prefixed after the add -- it does not need to be zeroed */ 1294 #if DOUBLE 1295 #define COFF 4 /* offset into acc */ 1296 #elif QUAD 1297 UBFROMUS(acc+4, 0); /* prefix 00 */ 1298 #define COFF 6 /* offset into acc */ 1299 #endif 1300 1301 GETCOEFF(dfl, acc+COFF); /* decode from decFloat */ 1302 ulsd=acc+COFF+DECPMAX-1; 1303 umsd=acc+4; /* [having this here avoids */ 1304 1305 #if DECTRACE 1306 {bcdnum tum; 1307 tum.msd=umsd; 1308 tum.lsd=ulsd; 1309 tum.exponent=bexpl-DECBIAS; 1310 tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign; 1311 decShowNum(&tum, "dflx");} 1312 #endif 1313 1314 /* if signs differ, take ten's complement of lhs (here the */ 1315 /* coefficient is subtracted from all-nines; the 1 is added during */ 1316 /* the later add cycle -- zeros to the right do not matter because */ 1317 /* the complement of zero is zero); these are fixed-length inverts */ 1318 /* where the lsd is known to be at a 4-byte boundary (so no borrow */ 1319 /* possible) */ 1320 carry=0; /* assume no carry */ 1321 if (diffsign) { 1322 carry=CARRYPAT; /* for +1 during add */ 1323 UBFROMUI(acc+ 4, 0x09090909-UBTOUI(acc+ 4)); 1324 UBFROMUI(acc+ 8, 0x09090909-UBTOUI(acc+ 8)); 1325 UBFROMUI(acc+12, 0x09090909-UBTOUI(acc+12)); 1326 UBFROMUI(acc+16, 0x09090909-UBTOUI(acc+16)); 1327 #if QUAD 1328 UBFROMUI(acc+20, 0x09090909-UBTOUI(acc+20)); 1329 UBFROMUI(acc+24, 0x09090909-UBTOUI(acc+24)); 1330 UBFROMUI(acc+28, 0x09090909-UBTOUI(acc+28)); 1331 UBFROMUI(acc+32, 0x09090909-UBTOUI(acc+32)); 1332 UBFROMUI(acc+36, 0x09090909-UBTOUI(acc+36)); 1333 #endif 1334 } /* diffsign */ 1335 1336 /* now process the rhs coefficient; if it cannot overlap lhs then */ 1337 /* it can be put straight into acc (with an appropriate gap, if */ 1338 /* needed) because no actual addition will be needed (except */ 1339 /* possibly to complete ten's complement) */ 1340 overlap=DECPMAX-(bexpl-bexpr); 1341 #if DECTRACE 1342 printf("exps: %ld %ld\n", (LI)(bexpl-DECBIAS), (LI)(bexpr-DECBIAS)); 1343 printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry); 1344 #endif 1345 1346 if (overlap<=0) { /* no overlap possible */ 1347 uInt gap; /* local work */ 1348 /* since a full addition is not needed, a ten's complement */ 1349 /* calculation started above may need to be completed */ 1350 if (carry) { 1351 for (ub=ulsd; *ub==9; ub--) *ub=0; 1352 *ub+=1; 1353 carry=0; /* taken care of */ 1354 } 1355 /* up to DECPMAX-1 digits of the final result can extend down */ 1356 /* below the LSD of the lhs, so if the gap is >DECPMAX then the */ 1357 /* rhs will be simply sticky bits. In this case the gap is */ 1358 /* clamped to DECPMAX and the exponent adjusted to suit [this is */ 1359 /* safe because the lhs is non-zero]. */ 1360 gap=-overlap; 1361 if (gap>DECPMAX) { 1362 bexpr+=gap-1; 1363 gap=DECPMAX; 1364 } 1365 ub=ulsd+gap+1; /* where MSD will go */ 1366 /* Fill the gap with 0s; note that there is no addition to do */ 1367 ut=acc+COFF+DECPMAX; /* start of gap */ 1368 for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* mind the gap */ 1369 if (overlap<-DECPMAX) { /* gap was > DECPMAX */ 1370 *ub=(uByte)(!DFISZERO(dfr)); /* make sticky digit */ 1371 } 1372 else { /* need full coefficient */ 1373 GETCOEFF(dfr, ub); /* decode from decFloat */ 1374 ub+=DECPMAX-1; /* new LSD... */ 1375 } 1376 ulsd=ub; /* save new LSD */ 1377 } /* no overlap possible */ 1378 1379 else { /* overlap>0 */ 1380 /* coefficients overlap (perhaps completely, although also */ 1381 /* perhaps only where zeros) */ 1382 if (overlap==DECPMAX) { /* aligned */ 1383 ub=buf+COFF; /* where msd will go */ 1384 #if QUAD 1385 UBFROMUS(buf+4, 0); /* clear quad's 00 */ 1386 #endif 1387 GETCOEFF(dfr, ub); /* decode from decFloat */ 1388 } 1389 else { /* unaligned */ 1390 ub=buf+COFF+DECPMAX-overlap; /* where MSD will go */ 1391 /* Fill the prefix gap with 0s; 8 will cover most common */ 1392 /* unalignments, so start with direct assignments (a loop is */ 1393 /* then used for any remaining -- the loop (and the one in a */ 1394 /* moment) is not then on the critical path because the number */ 1395 /* of additions is reduced by (at least) two in this case) */ 1396 UBFROMUI(buf+4, 0); /* [clears decQuad 00 too] */ 1397 UBFROMUI(buf+8, 0); 1398 if (ub>buf+12) { 1399 ut=buf+12; /* start any remaining */ 1400 for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* fill them */ 1401 } 1402 GETCOEFF(dfr, ub); /* decode from decFloat */ 1403 1404 /* now move tail of rhs across to main acc; again use direct */ 1405 /* copies for 8 digits-worth */ 1406 UBFROMUI(acc+COFF+DECPMAX, UBTOUI(buf+COFF+DECPMAX)); 1407 UBFROMUI(acc+COFF+DECPMAX+4, UBTOUI(buf+COFF+DECPMAX+4)); 1408 if (buf+COFF+DECPMAX+8<ub+DECPMAX) { 1409 us=buf+COFF+DECPMAX+8; /* source */ 1410 ut=acc+COFF+DECPMAX+8; /* target */ 1411 for (; us<ub+DECPMAX; us+=4, ut+=4) UBFROMUI(ut, UBTOUI(us)); 1412 } 1413 } /* unaligned */ 1414 1415 ulsd=acc+(ub-buf+DECPMAX-1); /* update LSD pointer */ 1416 1417 /* Now do the add of the non-tail; this is all nicely aligned, */ 1418 /* and is over a multiple of four digits (because for Quad two */ 1419 /* zero digits were added on the left); words in both acc and */ 1420 /* buf (buf especially) will often be zero */ 1421 /* [byte-by-byte add, here, is about 15% slower total effect than */ 1422 /* the by-fours] */ 1423 1424 /* Now effect the add; this is harder on a little-endian */ 1425 /* machine as the inter-digit carry cannot use the usual BCD */ 1426 /* addition trick because the bytes are loaded in the wrong order */ 1427 /* [this loop could be unrolled, but probably scarcely worth it] */ 1428 1429 ut=acc+COFF+DECPMAX-4; /* target LSW (acc) */ 1430 us=buf+COFF+DECPMAX-4; /* source LSW (buf, to add to acc) */ 1431 1432 #if !DECLITEND 1433 for (; ut>=acc+4; ut-=4, us-=4) { /* big-endian add loop */ 1434 /* bcd8 add */ 1435 carry+=UBTOUI(us); /* rhs + carry */ 1436 if (carry==0) continue; /* no-op */ 1437 carry+=UBTOUI(ut); /* lhs */ 1438 /* Big-endian BCD adjust (uses internal carry) */ 1439 carry+=0x76f6f6f6; /* note top nibble not all bits */ 1440 /* apply BCD adjust and save */ 1441 UBFROMUI(ut, (carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4)); 1442 carry>>=31; /* true carry was at far left */ 1443 } /* add loop */ 1444 #else 1445 for (; ut>=acc+4; ut-=4, us-=4) { /* little-endian add loop */ 1446 /* bcd8 add */ 1447 carry+=UBTOUI(us); /* rhs + carry */ 1448 if (carry==0) continue; /* no-op [common if unaligned] */ 1449 carry+=UBTOUI(ut); /* lhs */ 1450 /* Little-endian BCD adjust; inter-digit carry must be manual */ 1451 /* because the lsb from the array will be in the most-significant */ 1452 /* byte of carry */ 1453 carry+=0x76767676; /* note no inter-byte carries */ 1454 carry+=(carry & 0x80000000)>>15; 1455 carry+=(carry & 0x00800000)>>15; 1456 carry+=(carry & 0x00008000)>>15; 1457 carry-=(carry & 0x60606060)>>4; /* BCD adjust back */ 1458 UBFROMUI(ut, carry & 0x0f0f0f0f); /* clear debris and save */ 1459 /* here, final carry-out bit is at 0x00000080; move it ready */ 1460 /* for next word-add (i.e., to 0x01000000) */ 1461 carry=(carry & 0x00000080)<<17; 1462 } /* add loop */ 1463 #endif 1464 1465 #if DECTRACE 1466 {bcdnum tum; 1467 printf("Add done, carry=%08lx, diffsign=%ld\n", (LI)carry, (LI)diffsign); 1468 tum.msd=umsd; /* acc+4; */ 1469 tum.lsd=ulsd; 1470 tum.exponent=0; 1471 tum.sign=0; 1472 decShowNum(&tum, "dfadd");} 1473 #endif 1474 } /* overlap possible */ 1475 1476 /* ordering here is a little strange in order to have slowest path */ 1477 /* first in GCC asm listing */ 1478 if (diffsign) { /* subtraction */ 1479 if (!carry) { /* no carry out means RHS<LHS */ 1480 /* borrowed -- take ten's complement */ 1481 /* sign is lhs sign */ 1482 num.sign=DFWORD(dfl, 0) & DECFLOAT_Sign; 1483 1484 /* invert the coefficient first by fours, then add one; space */ 1485 /* at the end of the buffer ensures the by-fours is always */ 1486 /* safe, but lsd+1 must be cleared to prevent a borrow */ 1487 /* if big-endian */ 1488 #if !DECLITEND 1489 *(ulsd+1)=0; 1490 #endif 1491 /* there are always at least four coefficient words */ 1492 UBFROMUI(umsd, 0x09090909-UBTOUI(umsd)); 1493 UBFROMUI(umsd+4, 0x09090909-UBTOUI(umsd+4)); 1494 UBFROMUI(umsd+8, 0x09090909-UBTOUI(umsd+8)); 1495 UBFROMUI(umsd+12, 0x09090909-UBTOUI(umsd+12)); 1496 #if DOUBLE 1497 #define BNEXT 16 1498 #elif QUAD 1499 UBFROMUI(umsd+16, 0x09090909-UBTOUI(umsd+16)); 1500 UBFROMUI(umsd+20, 0x09090909-UBTOUI(umsd+20)); 1501 UBFROMUI(umsd+24, 0x09090909-UBTOUI(umsd+24)); 1502 UBFROMUI(umsd+28, 0x09090909-UBTOUI(umsd+28)); 1503 UBFROMUI(umsd+32, 0x09090909-UBTOUI(umsd+32)); 1504 #define BNEXT 36 1505 #endif 1506 if (ulsd>=umsd+BNEXT) { /* unaligned */ 1507 /* eight will handle most unaligments for Double; 16 for Quad */ 1508 UBFROMUI(umsd+BNEXT, 0x09090909-UBTOUI(umsd+BNEXT)); 1509 UBFROMUI(umsd+BNEXT+4, 0x09090909-UBTOUI(umsd+BNEXT+4)); 1510 #if DOUBLE 1511 #define BNEXTY (BNEXT+8) 1512 #elif QUAD 1513 UBFROMUI(umsd+BNEXT+8, 0x09090909-UBTOUI(umsd+BNEXT+8)); 1514 UBFROMUI(umsd+BNEXT+12, 0x09090909-UBTOUI(umsd+BNEXT+12)); 1515 #define BNEXTY (BNEXT+16) 1516 #endif 1517 if (ulsd>=umsd+BNEXTY) { /* very unaligned */ 1518 ut=umsd+BNEXTY; /* -> continue */ 1519 for (;;ut+=4) { 1520 UBFROMUI(ut, 0x09090909-UBTOUI(ut)); /* invert four digits */ 1521 if (ut>=ulsd-3) break; /* all done */ 1522 } 1523 } 1524 } 1525 /* complete the ten's complement by adding 1 */ 1526 for (ub=ulsd; *ub==9; ub--) *ub=0; 1527 *ub+=1; 1528 } /* borrowed */ 1529 1530 else { /* carry out means RHS>=LHS */ 1531 num.sign=DFWORD(dfr, 0) & DECFLOAT_Sign; 1532 /* all done except for the special IEEE 754 exact-zero-result */ 1533 /* rule (see above); while testing for zero, strip leading */ 1534 /* zeros (which will save decFinalize doing it) (this is in */ 1535 /* diffsign path, so carry impossible and true umsd is */ 1536 /* acc+COFF) */ 1537 1538 /* Check the initial coefficient area using the fast macro; */ 1539 /* this will often be all that needs to be done (as on the */ 1540 /* worst-case path when the subtraction was aligned and */ 1541 /* full-length) */ 1542 if (ISCOEFFZERO(acc+COFF)) { 1543 umsd=acc+COFF+DECPMAX-1; /* so far, so zero */ 1544 if (ulsd>umsd) { /* more to check */ 1545 umsd++; /* to align after checked area */ 1546 for (; UBTOUI(umsd)==0 && umsd+3<ulsd;) umsd+=4; 1547 for (; *umsd==0 && umsd<ulsd;) umsd++; 1548 } 1549 if (*umsd==0) { /* must be true zero (and diffsign) */ 1550 num.sign=0; /* assume + */ 1551 if (set->round==DEC_ROUND_FLOOR) num.sign=DECFLOAT_Sign; 1552 } 1553 } 1554 /* [else was not zero, might still have leading zeros] */ 1555 } /* subtraction gave positive result */ 1556 } /* diffsign */ 1557 1558 else { /* same-sign addition */ 1559 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign; 1560 #if DOUBLE 1561 if (carry) { /* only possible with decDouble */ 1562 *(acc+3)=1; /* [Quad has leading 00] */ 1563 umsd=acc+3; 1564 } 1565 #endif 1566 } /* same sign */ 1567 1568 num.msd=umsd; /* set MSD .. */ 1569 num.lsd=ulsd; /* .. and LSD */ 1570 num.exponent=bexpr-DECBIAS; /* set exponent to smaller, unbiassed */ 1571 1572 #if DECTRACE 1573 decFloatShow(dfl, "dfl"); 1574 decFloatShow(dfr, "dfr"); 1575 decShowNum(&num, "postadd"); 1576 #endif 1577 return decFinalize(result, &num, set); /* round, check, and lay out */ 1578 } /* decFloatAdd */ 1579 1580/* ------------------------------------------------------------------ */ 1581/* decFloatAnd -- logical digitwise AND of two decFloats */ 1582/* */ 1583/* result gets the result of ANDing dfl and dfr */ 1584/* dfl is the first decFloat (lhs) */ 1585/* dfr is the second decFloat (rhs) */ 1586/* set is the context */ 1587/* returns result, which will be canonical with sign=0 */ 1588/* */ 1589/* The operands must be positive, finite with exponent q=0, and */ 1590/* comprise just zeros and ones; if not, Invalid operation results. */ 1591/* ------------------------------------------------------------------ */ 1592decFloat * decFloatAnd(decFloat *result, 1593 const decFloat *dfl, const decFloat *dfr, 1594 decContext *set) { 1595 if (!DFISUINT01(dfl) || !DFISUINT01(dfr) 1596 || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set); 1597 /* the operands are positive finite integers (q=0) with just 0s and 1s */ 1598 #if DOUBLE 1599 DFWORD(result, 0)=ZEROWORD 1600 |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04009124); 1601 DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x49124491; 1602 #elif QUAD 1603 DFWORD(result, 0)=ZEROWORD 1604 |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04000912); 1605 DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x44912449; 1606 DFWORD(result, 2)=(DFWORD(dfl, 2) & DFWORD(dfr, 2))&0x12449124; 1607 DFWORD(result, 3)=(DFWORD(dfl, 3) & DFWORD(dfr, 3))&0x49124491; 1608 #endif 1609 return result; 1610 } /* decFloatAnd */ 1611 1612/* ------------------------------------------------------------------ */ 1613/* decFloatCanonical -- copy a decFloat, making canonical */ 1614/* */ 1615/* result gets the canonicalized df */ 1616/* df is the decFloat to copy and make canonical */ 1617/* returns result */ 1618/* */ 1619/* This works on specials, too; no error or exception is possible. */ 1620/* ------------------------------------------------------------------ */ 1621decFloat * decFloatCanonical(decFloat *result, const decFloat *df) { 1622 return decCanonical(result, df); 1623 } /* decFloatCanonical */ 1624 1625/* ------------------------------------------------------------------ */ 1626/* decFloatClass -- return the class of a decFloat */ 1627/* */ 1628/* df is the decFloat to test */ 1629/* returns the decClass that df falls into */ 1630/* ------------------------------------------------------------------ */ 1631enum decClass decFloatClass(const decFloat *df) { 1632 Int exp; /* exponent */ 1633 if (DFISSPECIAL(df)) { 1634 if (DFISQNAN(df)) return DEC_CLASS_QNAN; 1635 if (DFISSNAN(df)) return DEC_CLASS_SNAN; 1636 /* must be an infinity */ 1637 if (DFISSIGNED(df)) return DEC_CLASS_NEG_INF; 1638 return DEC_CLASS_POS_INF; 1639 } 1640 if (DFISZERO(df)) { /* quite common */ 1641 if (DFISSIGNED(df)) return DEC_CLASS_NEG_ZERO; 1642 return DEC_CLASS_POS_ZERO; 1643 } 1644 /* is finite and non-zero; similar code to decFloatIsNormal, here */ 1645 /* [this could be speeded up slightly by in-lining decFloatDigits] */ 1646 exp=GETEXPUN(df) /* get unbiased exponent .. */ 1647 +decFloatDigits(df)-1; /* .. and make adjusted exponent */ 1648 if (exp>=DECEMIN) { /* is normal */ 1649 if (DFISSIGNED(df)) return DEC_CLASS_NEG_NORMAL; 1650 return DEC_CLASS_POS_NORMAL; 1651 } 1652 /* is subnormal */ 1653 if (DFISSIGNED(df)) return DEC_CLASS_NEG_SUBNORMAL; 1654 return DEC_CLASS_POS_SUBNORMAL; 1655 } /* decFloatClass */ 1656 1657/* ------------------------------------------------------------------ */ 1658/* decFloatClassString -- return the class of a decFloat as a string */ 1659/* */ 1660/* df is the decFloat to test */ 1661/* returns a constant string describing the class df falls into */ 1662/* ------------------------------------------------------------------ */ 1663const char *decFloatClassString(const decFloat *df) { 1664 enum decClass eclass=decFloatClass(df); 1665 if (eclass==DEC_CLASS_POS_NORMAL) return DEC_ClassString_PN; 1666 if (eclass==DEC_CLASS_NEG_NORMAL) return DEC_ClassString_NN; 1667 if (eclass==DEC_CLASS_POS_ZERO) return DEC_ClassString_PZ; 1668 if (eclass==DEC_CLASS_NEG_ZERO) return DEC_ClassString_NZ; 1669 if (eclass==DEC_CLASS_POS_SUBNORMAL) return DEC_ClassString_PS; 1670 if (eclass==DEC_CLASS_NEG_SUBNORMAL) return DEC_ClassString_NS; 1671 if (eclass==DEC_CLASS_POS_INF) return DEC_ClassString_PI; 1672 if (eclass==DEC_CLASS_NEG_INF) return DEC_ClassString_NI; 1673 if (eclass==DEC_CLASS_QNAN) return DEC_ClassString_QN; 1674 if (eclass==DEC_CLASS_SNAN) return DEC_ClassString_SN; 1675 return DEC_ClassString_UN; /* Unknown */ 1676 } /* decFloatClassString */ 1677 1678/* ------------------------------------------------------------------ */ 1679/* decFloatCompare -- compare two decFloats; quiet NaNs allowed */ 1680/* */ 1681/* result gets the result of comparing dfl and dfr */ 1682/* dfl is the first decFloat (lhs) */ 1683/* dfr is the second decFloat (rhs) */ 1684/* set is the context */ 1685/* returns result, which may be -1, 0, 1, or NaN (Unordered) */ 1686/* ------------------------------------------------------------------ */ 1687decFloat * decFloatCompare(decFloat *result, 1688 const decFloat *dfl, const decFloat *dfr, 1689 decContext *set) { 1690 Int comp; /* work */ 1691 /* NaNs are handled as usual */ 1692 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set); 1693 /* numeric comparison needed */ 1694 comp=decNumCompare(dfl, dfr, 0); 1695 decFloatZero(result); 1696 if (comp==0) return result; 1697 DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */ 1698 if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */ 1699 return result; 1700 } /* decFloatCompare */ 1701 1702/* ------------------------------------------------------------------ */ 1703/* decFloatCompareSignal -- compare two decFloats; all NaNs signal */ 1704/* */ 1705/* result gets the result of comparing dfl and dfr */ 1706/* dfl is the first decFloat (lhs) */ 1707/* dfr is the second decFloat (rhs) */ 1708/* set is the context */ 1709/* returns result, which may be -1, 0, 1, or NaN (Unordered) */ 1710/* ------------------------------------------------------------------ */ 1711decFloat * decFloatCompareSignal(decFloat *result, 1712 const decFloat *dfl, const decFloat *dfr, 1713 decContext *set) { 1714 Int comp; /* work */ 1715 /* NaNs are handled as usual, except that all NaNs signal */ 1716 if (DFISNAN(dfl) || DFISNAN(dfr)) { 1717 set->status|=DEC_Invalid_operation; 1718 return decNaNs(result, dfl, dfr, set); 1719 } 1720 /* numeric comparison needed */ 1721 comp=decNumCompare(dfl, dfr, 0); 1722 decFloatZero(result); 1723 if (comp==0) return result; 1724 DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */ 1725 if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */ 1726 return result; 1727 } /* decFloatCompareSignal */ 1728 1729/* ------------------------------------------------------------------ */ 1730/* decFloatCompareTotal -- compare two decFloats with total ordering */ 1731/* */ 1732/* result gets the result of comparing dfl and dfr */ 1733/* dfl is the first decFloat (lhs) */ 1734/* dfr is the second decFloat (rhs) */ 1735/* returns result, which may be -1, 0, or 1 */ 1736/* ------------------------------------------------------------------ */ 1737decFloat * decFloatCompareTotal(decFloat *result, 1738 const decFloat *dfl, const decFloat *dfr) { 1739 Int comp; /* work */ 1740 uInt uiwork; /* for macros */ 1741 #if QUAD 1742 uShort uswork; /* .. */ 1743 #endif 1744 if (DFISNAN(dfl) || DFISNAN(dfr)) { 1745 Int nanl, nanr; /* work */ 1746 /* morph NaNs to +/- 1 or 2, leave numbers as 0 */ 1747 nanl=DFISSNAN(dfl)+DFISQNAN(dfl)*2; /* quiet > signalling */ 1748 if (DFISSIGNED(dfl)) nanl=-nanl; 1749 nanr=DFISSNAN(dfr)+DFISQNAN(dfr)*2; 1750 if (DFISSIGNED(dfr)) nanr=-nanr; 1751 if (nanl>nanr) comp=+1; 1752 else if (nanl<nanr) comp=-1; 1753 else { /* NaNs are the same type and sign .. must compare payload */ 1754 /* buffers need +2 for QUAD */ 1755 uByte bufl[DECPMAX+4]; /* for LHS coefficient + foot */ 1756 uByte bufr[DECPMAX+4]; /* for RHS coefficient + foot */ 1757 uByte *ub, *uc; /* work */ 1758 Int sigl; /* signum of LHS */ 1759 sigl=(DFISSIGNED(dfl) ? -1 : +1); 1760 1761 /* decode the coefficients */ 1762 /* (shift both right two if Quad to make a multiple of four) */ 1763 #if QUAD 1764 UBFROMUS(bufl, 0); 1765 UBFROMUS(bufr, 0); 1766 #endif 1767 GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */ 1768 GETCOEFF(dfr, bufr+QUAD*2); /* .. */ 1769 /* all multiples of four, here */ 1770 comp=0; /* assume equal */ 1771 for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) { 1772 uInt ui=UBTOUI(ub); 1773 if (ui==UBTOUI(uc)) continue; /* so far so same */ 1774 /* about to find a winner; go by bytes in case little-endian */ 1775 for (;; ub++, uc++) { 1776 if (*ub==*uc) continue; 1777 if (*ub>*uc) comp=sigl; /* difference found */ 1778 else comp=-sigl; /* .. */ 1779 break; 1780 } 1781 } 1782 } /* same NaN type and sign */ 1783 } 1784 else { 1785 /* numeric comparison needed */ 1786 comp=decNumCompare(dfl, dfr, 1); /* total ordering */ 1787 } 1788 decFloatZero(result); 1789 if (comp==0) return result; 1790 DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */ 1791 if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */ 1792 return result; 1793 } /* decFloatCompareTotal */ 1794 1795/* ------------------------------------------------------------------ */ 1796/* decFloatCompareTotalMag -- compare magnitudes with total ordering */ 1797/* */ 1798/* result gets the result of comparing abs(dfl) and abs(dfr) */ 1799/* dfl is the first decFloat (lhs) */ 1800/* dfr is the second decFloat (rhs) */ 1801/* returns result, which may be -1, 0, or 1 */ 1802/* ------------------------------------------------------------------ */ 1803decFloat * decFloatCompareTotalMag(decFloat *result, 1804 const decFloat *dfl, const decFloat *dfr) { 1805 decFloat a, b; /* for copy if needed */ 1806 /* copy and redirect signed operand(s) */ 1807 if (DFISSIGNED(dfl)) { 1808 decFloatCopyAbs(&a, dfl); 1809 dfl=&a; 1810 } 1811 if (DFISSIGNED(dfr)) { 1812 decFloatCopyAbs(&b, dfr); 1813 dfr=&b; 1814 } 1815 return decFloatCompareTotal(result, dfl, dfr); 1816 } /* decFloatCompareTotalMag */ 1817 1818/* ------------------------------------------------------------------ */ 1819/* decFloatCopy -- copy a decFloat as-is */ 1820/* */ 1821/* result gets the copy of dfl */ 1822/* dfl is the decFloat to copy */ 1823/* returns result */ 1824/* */ 1825/* This is a bitwise operation; no errors or exceptions are possible. */ 1826/* ------------------------------------------------------------------ */ 1827decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) { 1828 if (dfl!=result) *result=*dfl; /* copy needed */ 1829 return result; 1830 } /* decFloatCopy */ 1831 1832/* ------------------------------------------------------------------ */ 1833/* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0 */ 1834/* */ 1835/* result gets the copy of dfl with sign bit 0 */ 1836/* dfl is the decFloat to copy */ 1837/* returns result */ 1838/* */ 1839/* This is a bitwise operation; no errors or exceptions are possible. */ 1840/* ------------------------------------------------------------------ */ 1841decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) { 1842 if (dfl!=result) *result=*dfl; /* copy needed */ 1843 DFBYTE(result, 0)&=~0x80; /* zero sign bit */ 1844 return result; 1845 } /* decFloatCopyAbs */ 1846 1847/* ------------------------------------------------------------------ */ 1848/* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */ 1849/* */ 1850/* result gets the copy of dfl with sign bit inverted */ 1851/* dfl is the decFloat to copy */ 1852/* returns result */ 1853/* */ 1854/* This is a bitwise operation; no errors or exceptions are possible. */ 1855/* ------------------------------------------------------------------ */ 1856decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) { 1857 if (dfl!=result) *result=*dfl; /* copy needed */ 1858 DFBYTE(result, 0)^=0x80; /* invert sign bit */ 1859 return result; 1860 } /* decFloatCopyNegate */ 1861 1862/* ------------------------------------------------------------------ */ 1863/* decFloatCopySign -- copy a decFloat with the sign of another */ 1864/* */ 1865/* result gets the result of copying dfl with the sign of dfr */ 1866/* dfl is the first decFloat (lhs) */ 1867/* dfr is the second decFloat (rhs) */ 1868/* returns result */ 1869/* */ 1870/* This is a bitwise operation; no errors or exceptions are possible. */ 1871/* ------------------------------------------------------------------ */ 1872decFloat * decFloatCopySign(decFloat *result, 1873 const decFloat *dfl, const decFloat *dfr) { 1874 uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80); /* save sign bit */ 1875 if (dfl!=result) *result=*dfl; /* copy needed */ 1876 DFBYTE(result, 0)&=~0x80; /* clear sign .. */ 1877 DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* .. and set saved */ 1878 return result; 1879 } /* decFloatCopySign */ 1880 1881/* ------------------------------------------------------------------ */ 1882/* decFloatDigits -- return the number of digits in a decFloat */ 1883/* */ 1884/* df is the decFloat to investigate */ 1885/* returns the number of significant digits in the decFloat; a */ 1886/* zero coefficient returns 1 as does an infinity (a NaN returns */ 1887/* the number of digits in the payload) */ 1888/* ------------------------------------------------------------------ */ 1889/* private macro to extract a declet according to provided formula */ 1890/* (form), and if it is non-zero then return the calculated digits */ 1891/* depending on the declet number (n), where n=0 for the most */ 1892/* significant declet; uses uInt dpd for work */ 1893#define dpdlenchk(n, form) {dpd=(form)&0x3ff; \ 1894 if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);} 1895/* next one is used when it is known that the declet must be */ 1896/* non-zero, or is the final zero declet */ 1897#define dpdlendun(n, form) {dpd=(form)&0x3ff; \ 1898 if (dpd==0) return 1; \ 1899 return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);} 1900 1901uInt decFloatDigits(const decFloat *df) { 1902 uInt dpd; /* work */ 1903 uInt sourhi=DFWORD(df, 0); /* top word from source decFloat */ 1904 #if QUAD 1905 uInt sourmh, sourml; 1906 #endif 1907 uInt sourlo; 1908 1909 if (DFISINF(df)) return 1; 1910 /* A NaN effectively has an MSD of 0; otherwise if non-zero MSD */ 1911 /* then the coefficient is full-length */ 1912 if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX; 1913 1914 #if DOUBLE 1915 if (sourhi&0x0003ffff) { /* ends in first */ 1916 dpdlenchk(0, sourhi>>8); 1917 sourlo=DFWORD(df, 1); 1918 dpdlendun(1, (sourhi<<2) | (sourlo>>30)); 1919 } /* [cannot drop through] */ 1920 sourlo=DFWORD(df, 1); /* sourhi not involved now */ 1921 if (sourlo&0xfff00000) { /* in one of first two */ 1922 dpdlenchk(1, sourlo>>30); /* very rare */ 1923 dpdlendun(2, sourlo>>20); 1924 } /* [cannot drop through] */ 1925 dpdlenchk(3, sourlo>>10); 1926 dpdlendun(4, sourlo); 1927 /* [cannot drop through] */ 1928 1929 #elif QUAD 1930 if (sourhi&0x00003fff) { /* ends in first */ 1931 dpdlenchk(0, sourhi>>4); 1932 sourmh=DFWORD(df, 1); 1933 dpdlendun(1, ((sourhi)<<6) | (sourmh>>26)); 1934 } /* [cannot drop through] */ 1935 sourmh=DFWORD(df, 1); 1936 if (sourmh) { 1937 dpdlenchk(1, sourmh>>26); 1938 dpdlenchk(2, sourmh>>16); 1939 dpdlenchk(3, sourmh>>6); 1940 sourml=DFWORD(df, 2); 1941 dpdlendun(4, ((sourmh)<<4) | (sourml>>28)); 1942 } /* [cannot drop through] */ 1943 sourml=DFWORD(df, 2); 1944 if (sourml) { 1945 dpdlenchk(4, sourml>>28); 1946 dpdlenchk(5, sourml>>18); 1947 dpdlenchk(6, sourml>>8); 1948 sourlo=DFWORD(df, 3); 1949 dpdlendun(7, ((sourml)<<2) | (sourlo>>30)); 1950 } /* [cannot drop through] */ 1951 sourlo=DFWORD(df, 3); 1952 if (sourlo&0xfff00000) { /* in one of first two */ 1953 dpdlenchk(7, sourlo>>30); /* very rare */ 1954 dpdlendun(8, sourlo>>20); 1955 } /* [cannot drop through] */ 1956 dpdlenchk(9, sourlo>>10); 1957 dpdlendun(10, sourlo); 1958 /* [cannot drop through] */ 1959 #endif 1960 } /* decFloatDigits */ 1961 1962/* ------------------------------------------------------------------ */ 1963/* decFloatDivide -- divide a decFloat by another */ 1964/* */ 1965/* result gets the result of dividing dfl by dfr: */ 1966/* dfl is the first decFloat (lhs) */ 1967/* dfr is the second decFloat (rhs) */ 1968/* set is the context */ 1969/* returns result */ 1970/* */ 1971/* ------------------------------------------------------------------ */ 1972/* This is just a wrapper. */ 1973decFloat * decFloatDivide(decFloat *result, 1974 const decFloat *dfl, const decFloat *dfr, 1975 decContext *set) { 1976 return decDivide(result, dfl, dfr, set, DIVIDE); 1977 } /* decFloatDivide */ 1978 1979/* ------------------------------------------------------------------ */ 1980/* decFloatDivideInteger -- integer divide a decFloat by another */ 1981/* */ 1982/* result gets the result of dividing dfl by dfr: */ 1983/* dfl is the first decFloat (lhs) */ 1984/* dfr is the second decFloat (rhs) */ 1985/* set is the context */ 1986/* returns result */ 1987/* */ 1988/* ------------------------------------------------------------------ */ 1989decFloat * decFloatDivideInteger(decFloat *result, 1990 const decFloat *dfl, const decFloat *dfr, 1991 decContext *set) { 1992 return decDivide(result, dfl, dfr, set, DIVIDEINT); 1993 } /* decFloatDivideInteger */ 1994 1995/* ------------------------------------------------------------------ */ 1996/* decFloatFMA -- multiply and add three decFloats, fused */ 1997/* */ 1998/* result gets the result of (dfl*dfr)+dff with a single rounding */ 1999/* dfl is the first decFloat (lhs) */ 2000/* dfr is the second decFloat (rhs) */ 2001/* dff is the final decFloat (fhs) */ 2002/* set is the context */ 2003/* returns result */ 2004/* */ 2005/* ------------------------------------------------------------------ */ 2006decFloat * decFloatFMA(decFloat *result, const decFloat *dfl, 2007 const decFloat *dfr, const decFloat *dff, 2008 decContext *set) { 2009 2010 /* The accumulator has the bytes needed for FiniteMultiply, plus */ 2011 /* one byte to the left in case of carry, plus DECPMAX+2 to the */ 2012 /* right for the final addition (up to full fhs + round & sticky) */ 2013 #define FMALEN (ROUNDUP4(1+ (DECPMAX9*18+1) +DECPMAX+2)) 2014 uByte acc[FMALEN]; /* for multiplied coefficient in BCD */ 2015 /* .. and for final result */ 2016 bcdnum mul; /* for multiplication result */ 2017 bcdnum fin; /* for final operand, expanded */ 2018 uByte coe[ROUNDUP4(DECPMAX)]; /* dff coefficient in BCD */ 2019 bcdnum *hi, *lo; /* bcdnum with higher/lower exponent */ 2020 uInt diffsign; /* non-zero if signs differ */ 2021 uInt hipad; /* pad digit for hi if needed */ 2022 Int padding; /* excess exponent */ 2023 uInt carry; /* +1 for ten's complement and during add */ 2024 uByte *ub, *uh, *ul; /* work */ 2025 uInt uiwork; /* for macros */ 2026 2027 /* handle all the special values [any special operand leads to a */ 2028 /* special result] */ 2029 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr) || DFISSPECIAL(dff)) { 2030 decFloat proxy; /* multiplication result proxy */ 2031 /* NaNs are handled as usual, giving priority to sNaNs */ 2032 if (DFISSNAN(dfl) || DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set); 2033 if (DFISSNAN(dff)) return decNaNs(result, dff, NULL, set); 2034 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set); 2035 if (DFISNAN(dff)) return decNaNs(result, dff, NULL, set); 2036 /* One or more of the three is infinite */ 2037 /* infinity times zero is bad */ 2038 decFloatZero(&proxy); 2039 if (DFISINF(dfl)) { 2040 if (DFISZERO(dfr)) return decInvalid(result, set); 2041 decInfinity(&proxy, &proxy); 2042 } 2043 else if (DFISINF(dfr)) { 2044 if (DFISZERO(dfl)) return decInvalid(result, set); 2045 decInfinity(&proxy, &proxy); 2046 } 2047 /* compute sign of multiplication and place in proxy */ 2048 DFWORD(&proxy, 0)|=(DFWORD(dfl, 0)^DFWORD(dfr, 0))&DECFLOAT_Sign; 2049 if (!DFISINF(dff)) return decFloatCopy(result, &proxy); 2050 /* dff is Infinite */ 2051 if (!DFISINF(&proxy)) return decInfinity(result, dff); 2052 /* both sides of addition are infinite; different sign is bad */ 2053 if ((DFWORD(dff, 0)&DECFLOAT_Sign)!=(DFWORD(&proxy, 0)&DECFLOAT_Sign)) 2054 return decInvalid(result, set); 2055 return decFloatCopy(result, &proxy); 2056 } 2057 2058 /* Here when all operands are finite */ 2059 2060 /* First multiply dfl*dfr */ 2061 decFiniteMultiply(&mul, acc+1, dfl, dfr); 2062 /* The multiply is complete, exact and unbounded, and described in */ 2063 /* mul with the coefficient held in acc[1...] */ 2064 2065 /* now add in dff; the algorithm is essentially the same as */ 2066 /* decFloatAdd, but the code is different because the code there */ 2067 /* is highly optimized for adding two numbers of the same size */ 2068 fin.exponent=GETEXPUN(dff); /* get dff exponent and sign */ 2069 fin.sign=DFWORD(dff, 0)&DECFLOAT_Sign; 2070 diffsign=mul.sign^fin.sign; /* note if signs differ */ 2071 fin.msd=coe; 2072 fin.lsd=coe+DECPMAX-1; 2073 GETCOEFF(dff, coe); /* extract the coefficient */ 2074 2075 /* now set hi and lo so that hi points to whichever of mul and fin */ 2076 /* has the higher exponent and lo points to the other [don't care, */ 2077 /* if the same]. One coefficient will be in acc, the other in coe. */ 2078 if (mul.exponent>=fin.exponent) { 2079 hi=&mul; 2080 lo=&fin; 2081 } 2082 else { 2083 hi=&fin; 2084 lo=&mul; 2085 } 2086 2087 /* remove leading zeros on both operands; this will save time later */ 2088 /* and make testing for zero trivial (tests are safe because acc */ 2089 /* and coe are rounded up to uInts) */ 2090 for (; UBTOUI(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4; 2091 for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++; 2092 for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4; 2093 for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++; 2094 2095 /* if hi is zero then result will be lo (which has the smaller */ 2096 /* exponent), which also may need to be tested for zero for the */ 2097 /* weird IEEE 754 sign rules */ 2098 if (*hi->msd==0) { /* hi is zero */ 2099 /* "When the sum of two operands with opposite signs is */ 2100 /* exactly zero, the sign of that sum shall be '+' in all */ 2101 /* rounding modes except round toward -Infinity, in which */ 2102 /* mode that sign shall be '-'." */ 2103 if (diffsign) { 2104 if (*lo->msd==0) { /* lo is zero */ 2105 lo->sign=0; 2106 if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign; 2107 } /* diffsign && lo=0 */ 2108 } /* diffsign */ 2109 return decFinalize(result, lo, set); /* may need clamping */ 2110 } /* numfl is zero */ 2111 /* [here, both are minimal length and hi is non-zero] */ 2112 /* (if lo is zero then padding with zeros may be needed, below) */ 2113 2114 /* if signs differ, take the ten's complement of hi (zeros to the */ 2115 /* right do not matter because the complement of zero is zero); the */ 2116 /* +1 is done later, as part of the addition, inserted at the */ 2117 /* correct digit */ 2118 hipad=0; 2119 carry=0; 2120 if (diffsign) { 2121 hipad=9; 2122 carry=1; 2123 /* exactly the correct number of digits must be inverted */ 2124 for (uh=hi->msd; uh<hi->lsd-3; uh+=4) UBFROMUI(uh, 0x09090909-UBTOUI(uh)); 2125 for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh); 2126 } 2127 2128 /* ready to add; note that hi has no leading zeros so gap */ 2129 /* calculation does not have to be as pessimistic as in decFloatAdd */ 2130 /* (this is much more like the arbitrary-precision algorithm in */ 2131 /* Rexx and decNumber) */ 2132 2133 /* padding is the number of zeros that would need to be added to hi */ 2134 /* for its lsd to be aligned with the lsd of lo */ 2135 padding=hi->exponent-lo->exponent; 2136 /* printf("FMA pad %ld\n", (LI)padding); */ 2137 2138 /* the result of the addition will be built into the accumulator, */ 2139 /* starting from the far right; this could be either hi or lo, and */ 2140 /* will be aligned */ 2141 ub=acc+FMALEN-1; /* where lsd of result will go */ 2142 ul=lo->lsd; /* lsd of rhs */ 2143 2144 if (padding!=0) { /* unaligned */ 2145 /* if the msd of lo is more than DECPMAX+2 digits to the right of */ 2146 /* the original msd of hi then it can be reduced to a single */ 2147 /* digit at the right place, as it stays clear of hi digits */ 2148 /* [it must be DECPMAX+2 because during a subtraction the msd */ 2149 /* could become 0 after a borrow from 1.000 to 0.9999...] */ 2150 2151 Int hilen=(Int)(hi->lsd-hi->msd+1); /* length of hi */ 2152 Int lolen=(Int)(lo->lsd-lo->msd+1); /* and of lo */ 2153 2154 if (hilen+padding-lolen > DECPMAX+2) { /* can reduce lo to single */ 2155 /* make sure it is virtually at least DECPMAX from hi->msd, at */ 2156 /* least to right of hi->lsd (in case of destructive subtract), */ 2157 /* and separated by at least two digits from either of those */ 2158 /* (the tricky DOUBLE case is when hi is a 1 that will become a */ 2159 /* 0.9999... by subtraction: */ 2160 /* hi: 1 E+16 */ 2161 /* lo: .................1000000000000000 E-16 */ 2162 /* which for the addition pads to: */ 2163 /* hi: 1000000000000000000 E-16 */ 2164 /* lo: .................1000000000000000 E-16 */ 2165 Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3; 2166 2167 /* printf("FMA reduce: %ld\n", (LI)reduce); */ 2168 lo->lsd=lo->msd; /* to single digit [maybe 0] */ 2169 lo->exponent=newexp; /* new lowest exponent */ 2170 padding=hi->exponent-lo->exponent; /* recalculate */ 2171 ul=lo->lsd; /* .. and repoint */ 2172 } 2173 2174 /* padding is still > 0, but will fit in acc (less leading carry slot) */ 2175 #if DECCHECK 2176 if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding); 2177 if (hilen+padding+1>FMALEN) 2178 printf("FMA excess hilen+padding: %ld+%ld \n", (LI)hilen, (LI)padding); 2179 /* printf("FMA padding: %ld\n", (LI)padding); */ 2180 #endif 2181 2182 /* padding digits can now be set in the result; one or more of */ 2183 /* these will come from lo; others will be zeros in the gap */ 2184 for (; ul-3>=lo->msd && padding>3; padding-=4, ul-=4, ub-=4) { 2185 UBFROMUI(ub-3, UBTOUI(ul-3)); /* [cannot overlap] */ 2186 } 2187 for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul; 2188 for (;padding>0; padding--, ub--) *ub=0; /* mind the gap */ 2189 } 2190 2191 /* addition now complete to the right of the rightmost digit of hi */ 2192 uh=hi->lsd; 2193 2194 /* dow do the add from hi->lsd to the left */ 2195 /* [bytewise, because either operand can run out at any time] */ 2196 /* carry was set up depending on ten's complement above */ 2197 /* first assume both operands have some digits */ 2198 for (;; ub--) { 2199 if (uh<hi->msd || ul<lo->msd) break; 2200 *ub=(uByte)(carry+(*uh--)+(*ul--)); 2201 carry=0; 2202 if (*ub<10) continue; 2203 *ub-=10; 2204 carry=1; 2205 } /* both loop */ 2206 2207 if (ul<lo->msd) { /* to left of lo */ 2208 for (;; ub--) { 2209 if (uh<hi->msd) break; 2210 *ub=(uByte)(carry+(*uh--)); /* [+0] */ 2211 carry=0; 2212 if (*ub<10) continue; 2213 *ub-=10; 2214 carry=1; 2215 } /* hi loop */ 2216 } 2217 else { /* to left of hi */ 2218 for (;; ub--) { 2219 if (ul<lo->msd) break; 2220 *ub=(uByte)(carry+hipad+(*ul--)); 2221 carry=0; 2222 if (*ub<10) continue; 2223 *ub-=10; 2224 carry=1; 2225 } /* lo loop */ 2226 } 2227 2228 /* addition complete -- now handle carry, borrow, etc. */ 2229 /* use lo to set up the num (its exponent is already correct, and */ 2230 /* sign usually is) */ 2231 lo->msd=ub+1; 2232 lo->lsd=acc+FMALEN-1; 2233 /* decShowNum(lo, "lo"); */ 2234 if (!diffsign) { /* same-sign addition */ 2235 if (carry) { /* carry out */ 2236 *ub=1; /* place the 1 .. */ 2237 lo->msd--; /* .. and update */ 2238 } 2239 } /* same sign */ 2240 else { /* signs differed (subtraction) */ 2241 if (!carry) { /* no carry out means hi<lo */ 2242 /* borrowed -- take ten's complement of the right digits */ 2243 lo->sign=hi->sign; /* sign is lhs sign */ 2244 for (ul=lo->msd; ul<lo->lsd-3; ul+=4) UBFROMUI(ul, 0x09090909-UBTOUI(ul)); 2245 for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); /* [leaves ul at lsd+1] */ 2246 /* complete the ten's complement by adding 1 [cannot overrun] */ 2247 for (ul--; *ul==9; ul--) *ul=0; 2248 *ul+=1; 2249 } /* borrowed */ 2250 else { /* carry out means hi>=lo */ 2251 /* sign to use is lo->sign */ 2252 /* all done except for the special IEEE 754 exact-zero-result */ 2253 /* rule (see above); while testing for zero, strip leading */ 2254 /* zeros (which will save decFinalize doing it) */ 2255 for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4; 2256 for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++; 2257 if (*lo->msd==0) { /* must be true zero (and diffsign) */ 2258 lo->sign=0; /* assume + */ 2259 if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign; 2260 } 2261 /* [else was not zero, might still have leading zeros] */ 2262 } /* subtraction gave positive result */ 2263 } /* diffsign */ 2264 2265 #if DECCHECK 2266 /* assert no left underrun */ 2267 if (lo->msd<acc) { 2268 printf("FMA underrun by %ld \n", (LI)(acc-lo->msd)); 2269 } 2270 #endif 2271 2272 return decFinalize(result, lo, set); /* round, check, and lay out */ 2273 } /* decFloatFMA */ 2274 2275/* ------------------------------------------------------------------ */ 2276/* decFloatFromInt -- initialise a decFloat from an Int */ 2277/* */ 2278/* result gets the converted Int */ 2279/* n is the Int to convert */ 2280/* returns result */ 2281/* */ 2282/* The result is Exact; no errors or exceptions are possible. */ 2283/* ------------------------------------------------------------------ */ 2284decFloat * decFloatFromInt32(decFloat *result, Int n) { 2285 uInt u=(uInt)n; /* copy as bits */ 2286 uInt encode; /* work */ 2287 DFWORD(result, 0)=ZEROWORD; /* always */ 2288 #if QUAD 2289 DFWORD(result, 1)=0; 2290 DFWORD(result, 2)=0; 2291 #endif 2292 if (n<0) { /* handle -n with care */ 2293 /* [This can be done without the test, but is then slightly slower] */ 2294 u=(~u)+1; 2295 DFWORD(result, 0)|=DECFLOAT_Sign; 2296 } 2297 /* Since the maximum value of u now is 2**31, only the low word of */ 2298 /* result is affected */ 2299 encode=BIN2DPD[u%1000]; 2300 u/=1000; 2301 encode|=BIN2DPD[u%1000]<<10; 2302 u/=1000; 2303 encode|=BIN2DPD[u%1000]<<20; 2304 u/=1000; /* now 0, 1, or 2 */ 2305 encode|=u<<30; 2306 DFWORD(result, DECWORDS-1)=encode; 2307 return result; 2308 } /* decFloatFromInt32 */ 2309 2310/* ------------------------------------------------------------------ */ 2311/* decFloatFromUInt -- initialise a decFloat from a uInt */ 2312/* */ 2313/* result gets the converted uInt */ 2314/* n is the uInt to convert */ 2315/* returns result */ 2316/* */ 2317/* The result is Exact; no errors or exceptions are possible. */ 2318/* ------------------------------------------------------------------ */ 2319decFloat * decFloatFromUInt32(decFloat *result, uInt u) { 2320 uInt encode; /* work */ 2321 DFWORD(result, 0)=ZEROWORD; /* always */ 2322 #if QUAD 2323 DFWORD(result, 1)=0; 2324 DFWORD(result, 2)=0; 2325 #endif 2326 encode=BIN2DPD[u%1000]; 2327 u/=1000; 2328 encode|=BIN2DPD[u%1000]<<10; 2329 u/=1000; 2330 encode|=BIN2DPD[u%1000]<<20; 2331 u/=1000; /* now 0 -> 4 */ 2332 encode|=u<<30; 2333 DFWORD(result, DECWORDS-1)=encode; 2334 DFWORD(result, DECWORDS-2)|=u>>2; /* rarely non-zero */ 2335 return result; 2336 } /* decFloatFromUInt32 */ 2337 2338/* ------------------------------------------------------------------ */ 2339/* decFloatInvert -- logical digitwise INVERT of a decFloat */ 2340/* */ 2341/* result gets the result of INVERTing df */ 2342/* df is the decFloat to invert */ 2343/* set is the context */ 2344/* returns result, which will be canonical with sign=0 */ 2345/* */ 2346/* The operand must be positive, finite with exponent q=0, and */ 2347/* comprise just zeros and ones; if not, Invalid operation results. */ 2348/* ------------------------------------------------------------------ */ 2349decFloat * decFloatInvert(decFloat *result, const decFloat *df, 2350 decContext *set) { 2351 uInt sourhi=DFWORD(df, 0); /* top word of dfs */ 2352 2353 if (!DFISUINT01(df) || !DFISCC01(df)) return decInvalid(result, set); 2354 /* the operand is a finite integer (q=0) */ 2355 #if DOUBLE 2356 DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04009124); 2357 DFWORD(result, 1)=(~DFWORD(df, 1)) &0x49124491; 2358 #elif QUAD 2359 DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04000912); 2360 DFWORD(result, 1)=(~DFWORD(df, 1)) &0x44912449; 2361 DFWORD(result, 2)=(~DFWORD(df, 2)) &0x12449124; 2362 DFWORD(result, 3)=(~DFWORD(df, 3)) &0x49124491; 2363 #endif 2364 return result; 2365 } /* decFloatInvert */ 2366 2367/* ------------------------------------------------------------------ */ 2368/* decFloatIs -- decFloat tests (IsSigned, etc.) */ 2369/* */ 2370/* df is the decFloat to test */ 2371/* returns 0 or 1 in a uInt */ 2372/* */ 2373/* Many of these could be macros, but having them as real functions */ 2374/* is a little cleaner (and they can be referred to here by the */ 2375/* generic names) */ 2376/* ------------------------------------------------------------------ */ 2377uInt decFloatIsCanonical(const decFloat *df) { 2378 if (DFISSPECIAL(df)) { 2379 if (DFISINF(df)) { 2380 if (DFWORD(df, 0)&ECONMASK) return 0; /* exponent continuation */ 2381 if (!DFISCCZERO(df)) return 0; /* coefficient continuation */ 2382 return 1; 2383 } 2384 /* is a NaN */ 2385 if (DFWORD(df, 0)&ECONNANMASK) return 0; /* exponent continuation */ 2386 if (DFISCCZERO(df)) return 1; /* coefficient continuation */ 2387 /* drop through to check payload */ 2388 } 2389 { /* declare block */ 2390 #if DOUBLE 2391 uInt sourhi=DFWORD(df, 0); 2392 uInt sourlo=DFWORD(df, 1); 2393 if (CANONDPDOFF(sourhi, 8) 2394 && CANONDPDTWO(sourhi, sourlo, 30) 2395 && CANONDPDOFF(sourlo, 20) 2396 && CANONDPDOFF(sourlo, 10) 2397 && CANONDPDOFF(sourlo, 0)) return 1; 2398 #elif QUAD 2399 uInt sourhi=DFWORD(df, 0); 2400 uInt sourmh=DFWORD(df, 1); 2401 uInt sourml=DFWORD(df, 2); 2402 uInt sourlo=DFWORD(df, 3); 2403 if (CANONDPDOFF(sourhi, 4) 2404 && CANONDPDTWO(sourhi, sourmh, 26) 2405 && CANONDPDOFF(sourmh, 16) 2406 && CANONDPDOFF(sourmh, 6) 2407 && CANONDPDTWO(sourmh, sourml, 28) 2408 && CANONDPDOFF(sourml, 18) 2409 && CANONDPDOFF(sourml, 8) 2410 && CANONDPDTWO(sourml, sourlo, 30) 2411 && CANONDPDOFF(sourlo, 20) 2412 && CANONDPDOFF(sourlo, 10) 2413 && CANONDPDOFF(sourlo, 0)) return 1; 2414 #endif 2415 } /* block */ 2416 return 0; /* a declet is non-canonical */ 2417 } 2418 2419uInt decFloatIsFinite(const decFloat *df) { 2420 return !DFISSPECIAL(df); 2421 } 2422uInt decFloatIsInfinite(const decFloat *df) { 2423 return DFISINF(df); 2424 } 2425uInt decFloatIsInteger(const decFloat *df) { 2426 return DFISINT(df); 2427 } 2428uInt decFloatIsNaN(const decFloat *df) { 2429 return DFISNAN(df); 2430 } 2431uInt decFloatIsNormal(const decFloat *df) { 2432 Int exp; /* exponent */ 2433 if (DFISSPECIAL(df)) return 0; 2434 if (DFISZERO(df)) return 0; 2435 /* is finite and non-zero */ 2436 exp=GETEXPUN(df) /* get unbiased exponent .. */ 2437 +decFloatDigits(df)-1; /* .. and make adjusted exponent */ 2438 return (exp>=DECEMIN); /* < DECEMIN is subnormal */ 2439 } 2440uInt decFloatIsSignaling(const decFloat *df) { 2441 return DFISSNAN(df); 2442 } 2443uInt decFloatIsSignalling(const decFloat *df) { 2444 return DFISSNAN(df); 2445 } 2446uInt decFloatIsSigned(const decFloat *df) { 2447 return DFISSIGNED(df); 2448 } 2449uInt decFloatIsSubnormal(const decFloat *df) { 2450 if (DFISSPECIAL(df)) return 0; 2451 /* is finite */ 2452 if (decFloatIsNormal(df)) return 0; 2453 /* it is <Nmin, but could be zero */ 2454 if (DFISZERO(df)) return 0; 2455 return 1; /* is subnormal */ 2456 } 2457uInt decFloatIsZero(const decFloat *df) { 2458 return DFISZERO(df); 2459 } /* decFloatIs... */ 2460 2461/* ------------------------------------------------------------------ */ 2462/* decFloatLogB -- return adjusted exponent, by 754 rules */ 2463/* */ 2464/* result gets the adjusted exponent as an integer, or a NaN etc. */ 2465/* df is the decFloat to be examined */ 2466/* set is the context */ 2467/* returns result */ 2468/* */ 2469/* Notable cases: */ 2470/* A<0 -> Use |A| */ 2471/* A=0 -> -Infinity (Division by zero) */ 2472/* A=Infinite -> +Infinity (Exact) */ 2473/* A=1 exactly -> 0 (Exact) */ 2474/* NaNs are propagated as usual */ 2475/* ------------------------------------------------------------------ */ 2476decFloat * decFloatLogB(decFloat *result, const decFloat *df, 2477 decContext *set) { 2478 Int ae; /* adjusted exponent */ 2479 if (DFISNAN(df)) return decNaNs(result, df, NULL, set); 2480 if (DFISINF(df)) { 2481 DFWORD(result, 0)=0; /* need +ve */ 2482 return decInfinity(result, result); /* canonical +Infinity */ 2483 } 2484 if (DFISZERO(df)) { 2485 set->status|=DEC_Division_by_zero; /* as per 754 */ 2486 DFWORD(result, 0)=DECFLOAT_Sign; /* make negative */ 2487 return decInfinity(result, result); /* canonical -Infinity */ 2488 } 2489 ae=GETEXPUN(df) /* get unbiased exponent .. */ 2490 +decFloatDigits(df)-1; /* .. and make adjusted exponent */ 2491 /* ae has limited range (3 digits for DOUBLE and 4 for QUAD), so */ 2492 /* it is worth using a special case of decFloatFromInt32 */ 2493 DFWORD(result, 0)=ZEROWORD; /* always */ 2494 if (ae<0) { 2495 DFWORD(result, 0)|=DECFLOAT_Sign; /* -0 so far */ 2496 ae=-ae; 2497 } 2498 #if DOUBLE 2499 DFWORD(result, 1)=BIN2DPD[ae]; /* a single declet */ 2500 #elif QUAD 2501 DFWORD(result, 1)=0; 2502 DFWORD(result, 2)=0; 2503 DFWORD(result, 3)=(ae/1000)<<10; /* is <10, so need no DPD encode */ 2504 DFWORD(result, 3)|=BIN2DPD[ae%1000]; 2505 #endif 2506 return result; 2507 } /* decFloatLogB */ 2508 2509/* ------------------------------------------------------------------ */ 2510/* decFloatMax -- return maxnum of two operands */ 2511/* */ 2512/* result gets the chosen decFloat */ 2513/* dfl is the first decFloat (lhs) */ 2514/* dfr is the second decFloat (rhs) */ 2515/* set is the context */ 2516/* returns result */ 2517/* */ 2518/* If just one operand is a quiet NaN it is ignored. */ 2519/* ------------------------------------------------------------------ */ 2520decFloat * decFloatMax(decFloat *result, 2521 const decFloat *dfl, const decFloat *dfr, 2522 decContext *set) { 2523 Int comp; 2524 if (DFISNAN(dfl)) { 2525 /* sNaN or both NaNs leads to normal NaN processing */ 2526 if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set); 2527 return decCanonical(result, dfr); /* RHS is numeric */ 2528 } 2529 if (DFISNAN(dfr)) { 2530 /* sNaN leads to normal NaN processing (both NaN handled above) */ 2531 if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set); 2532 return decCanonical(result, dfl); /* LHS is numeric */ 2533 } 2534 /* Both operands are numeric; numeric comparison needed -- use */ 2535 /* total order for a well-defined choice (and +0 > -0) */ 2536 comp=decNumCompare(dfl, dfr, 1); 2537 if (comp>=0) return decCanonical(result, dfl); 2538 return decCanonical(result, dfr); 2539 } /* decFloatMax */ 2540 2541/* ------------------------------------------------------------------ */ 2542/* decFloatMaxMag -- return maxnummag of two operands */ 2543/* */ 2544/* result gets the chosen decFloat */ 2545/* dfl is the first decFloat (lhs) */ 2546/* dfr is the second decFloat (rhs) */ 2547/* set is the context */ 2548/* returns result */ 2549/* */ 2550/* Returns according to the magnitude comparisons if both numeric and */ 2551/* unequal, otherwise returns maxnum */ 2552/* ------------------------------------------------------------------ */ 2553decFloat * decFloatMaxMag(decFloat *result, 2554 const decFloat *dfl, const decFloat *dfr, 2555 decContext *set) { 2556 Int comp; 2557 decFloat absl, absr; 2558 if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMax(result, dfl, dfr, set); 2559 2560 decFloatCopyAbs(&absl, dfl); 2561 decFloatCopyAbs(&absr, dfr); 2562 comp=decNumCompare(&absl, &absr, 0); 2563 if (comp>0) return decCanonical(result, dfl); 2564 if (comp<0) return decCanonical(result, dfr); 2565 return decFloatMax(result, dfl, dfr, set); 2566 } /* decFloatMaxMag */ 2567 2568/* ------------------------------------------------------------------ */ 2569/* decFloatMin -- return minnum of two operands */ 2570/* */ 2571/* result gets the chosen decFloat */ 2572/* dfl is the first decFloat (lhs) */ 2573/* dfr is the second decFloat (rhs) */ 2574/* set is the context */ 2575/* returns result */ 2576/* */ 2577/* If just one operand is a quiet NaN it is ignored. */ 2578/* ------------------------------------------------------------------ */ 2579decFloat * decFloatMin(decFloat *result, 2580 const decFloat *dfl, const decFloat *dfr, 2581 decContext *set) { 2582 Int comp; 2583 if (DFISNAN(dfl)) { 2584 /* sNaN or both NaNs leads to normal NaN processing */ 2585 if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set); 2586 return decCanonical(result, dfr); /* RHS is numeric */ 2587 } 2588 if (DFISNAN(dfr)) { 2589 /* sNaN leads to normal NaN processing (both NaN handled above) */ 2590 if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set); 2591 return decCanonical(result, dfl); /* LHS is numeric */ 2592 } 2593 /* Both operands are numeric; numeric comparison needed -- use */ 2594 /* total order for a well-defined choice (and +0 > -0) */ 2595 comp=decNumCompare(dfl, dfr, 1); 2596 if (comp<=0) return decCanonical(result, dfl); 2597 return decCanonical(result, dfr); 2598 } /* decFloatMin */ 2599 2600/* ------------------------------------------------------------------ */ 2601/* decFloatMinMag -- return minnummag of two operands */ 2602/* */ 2603/* result gets the chosen decFloat */ 2604/* dfl is the first decFloat (lhs) */ 2605/* dfr is the second decFloat (rhs) */ 2606/* set is the context */ 2607/* returns result */ 2608/* */ 2609/* Returns according to the magnitude comparisons if both numeric and */ 2610/* unequal, otherwise returns minnum */ 2611/* ------------------------------------------------------------------ */ 2612decFloat * decFloatMinMag(decFloat *result, 2613 const decFloat *dfl, const decFloat *dfr, 2614 decContext *set) { 2615 Int comp; 2616 decFloat absl, absr; 2617 if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMin(result, dfl, dfr, set); 2618 2619 decFloatCopyAbs(&absl, dfl); 2620 decFloatCopyAbs(&absr, dfr); 2621 comp=decNumCompare(&absl, &absr, 0); 2622 if (comp<0) return decCanonical(result, dfl); 2623 if (comp>0) return decCanonical(result, dfr); 2624 return decFloatMin(result, dfl, dfr, set); 2625 } /* decFloatMinMag */ 2626 2627/* ------------------------------------------------------------------ */ 2628/* decFloatMinus -- negate value, heeding NaNs, etc. */ 2629/* */ 2630/* result gets the canonicalized 0-df */ 2631/* df is the decFloat to minus */ 2632/* set is the context */ 2633/* returns result */ 2634/* */ 2635/* This has the same effect as 0-df where the exponent of the zero is */ 2636/* the same as that of df (if df is finite). */ 2637/* The effect is also the same as decFloatCopyNegate except that NaNs */ 2638/* are handled normally (the sign of a NaN is not affected, and an */ 2639/* sNaN will signal), the result is canonical, and zero gets sign 0. */ 2640/* ------------------------------------------------------------------ */ 2641decFloat * decFloatMinus(decFloat *result, const decFloat *df, 2642 decContext *set) { 2643 if (DFISNAN(df)) return decNaNs(result, df, NULL, set); 2644 decCanonical(result, df); /* copy and check */ 2645 if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */ 2646 else DFBYTE(result, 0)^=0x80; /* flip sign bit */ 2647 return result; 2648 } /* decFloatMinus */ 2649 2650/* ------------------------------------------------------------------ */ 2651/* decFloatMultiply -- multiply two decFloats */ 2652/* */ 2653/* result gets the result of multiplying dfl and dfr: */ 2654/* dfl is the first decFloat (lhs) */ 2655/* dfr is the second decFloat (rhs) */ 2656/* set is the context */ 2657/* returns result */ 2658/* */ 2659/* ------------------------------------------------------------------ */ 2660decFloat * decFloatMultiply(decFloat *result, 2661 const decFloat *dfl, const decFloat *dfr, 2662 decContext *set) { 2663 bcdnum num; /* for final conversion */ 2664 uByte bcdacc[DECPMAX9*18+1]; /* for coefficent in BCD */ 2665 2666 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */ 2667 /* NaNs are handled as usual */ 2668 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set); 2669 /* infinity times zero is bad */ 2670 if (DFISINF(dfl) && DFISZERO(dfr)) return decInvalid(result, set); 2671 if (DFISINF(dfr) && DFISZERO(dfl)) return decInvalid(result, set); 2672 /* both infinite; return canonical infinity with computed sign */ 2673 DFWORD(result, 0)=DFWORD(dfl, 0)^DFWORD(dfr, 0); /* compute sign */ 2674 return decInfinity(result, result); 2675 } 2676 2677 /* Here when both operands are finite */ 2678 decFiniteMultiply(&num, bcdacc, dfl, dfr); 2679 return decFinalize(result, &num, set); /* round, check, and lay out */ 2680 } /* decFloatMultiply */ 2681 2682/* ------------------------------------------------------------------ */ 2683/* decFloatNextMinus -- next towards -Infinity */ 2684/* */ 2685/* result gets the next lesser decFloat */ 2686/* dfl is the decFloat to start with */ 2687/* set is the context */ 2688/* returns result */ 2689/* */ 2690/* This is 754 nextdown; Invalid is the only status possible (from */ 2691/* an sNaN). */ 2692/* ------------------------------------------------------------------ */ 2693decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl, 2694 decContext *set) { 2695 decFloat delta; /* tiny increment */ 2696 uInt savestat; /* saves status */ 2697 enum rounding saveround; /* .. and mode */ 2698 2699 /* +Infinity is the special case */ 2700 if (DFISINF(dfl) && !DFISSIGNED(dfl)) { 2701 DFSETNMAX(result); 2702 return result; /* [no status to set] */ 2703 } 2704 /* other cases are effected by sutracting a tiny delta -- this */ 2705 /* should be done in a wider format as the delta is unrepresentable */ 2706 /* here (but can be done with normal add if the sign of zero is */ 2707 /* treated carefully, because no Inexactitude is interesting); */ 2708 /* rounding to -Infinity then pushes the result to next below */ 2709 decFloatZero(&delta); /* set up tiny delta */ 2710 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */ 2711 DFWORD(&delta, 0)=DECFLOAT_Sign; /* Sign=1 + biased exponent=0 */ 2712 /* set up for the directional round */ 2713 saveround=set->round; /* save mode */ 2714 set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */ 2715 savestat=set->status; /* save status */ 2716 decFloatAdd(result, dfl, &delta, set); 2717 /* Add rules mess up the sign when going from +Ntiny to 0 */ 2718 if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */ 2719 set->status&=DEC_Invalid_operation; /* preserve only sNaN status */ 2720 set->status|=savestat; /* restore pending flags */ 2721 set->round=saveround; /* .. and mode */ 2722 return result; 2723 } /* decFloatNextMinus */ 2724 2725/* ------------------------------------------------------------------ */ 2726/* decFloatNextPlus -- next towards +Infinity */ 2727/* */ 2728/* result gets the next larger decFloat */ 2729/* dfl is the decFloat to start with */ 2730/* set is the context */ 2731/* returns result */ 2732/* */ 2733/* This is 754 nextup; Invalid is the only status possible (from */ 2734/* an sNaN). */ 2735/* ------------------------------------------------------------------ */ 2736decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl, 2737 decContext *set) { 2738 uInt savestat; /* saves status */ 2739 enum rounding saveround; /* .. and mode */ 2740 decFloat delta; /* tiny increment */ 2741 2742 /* -Infinity is the special case */ 2743 if (DFISINF(dfl) && DFISSIGNED(dfl)) { 2744 DFSETNMAX(result); 2745 DFWORD(result, 0)|=DECFLOAT_Sign; /* make negative */ 2746 return result; /* [no status to set] */ 2747 } 2748 /* other cases are effected by sutracting a tiny delta -- this */ 2749 /* should be done in a wider format as the delta is unrepresentable */ 2750 /* here (but can be done with normal add if the sign of zero is */ 2751 /* treated carefully, because no Inexactitude is interesting); */ 2752 /* rounding to +Infinity then pushes the result to next above */ 2753 decFloatZero(&delta); /* set up tiny delta */ 2754 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */ 2755 DFWORD(&delta, 0)=0; /* Sign=0 + biased exponent=0 */ 2756 /* set up for the directional round */ 2757 saveround=set->round; /* save mode */ 2758 set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */ 2759 savestat=set->status; /* save status */ 2760 decFloatAdd(result, dfl, &delta, set); 2761 /* Add rules mess up the sign when going from -Ntiny to -0 */ 2762 if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */ 2763 set->status&=DEC_Invalid_operation; /* preserve only sNaN status */ 2764 set->status|=savestat; /* restore pending flags */ 2765 set->round=saveround; /* .. and mode */ 2766 return result; 2767 } /* decFloatNextPlus */ 2768 2769/* ------------------------------------------------------------------ */ 2770/* decFloatNextToward -- next towards a decFloat */ 2771/* */ 2772/* result gets the next decFloat */ 2773/* dfl is the decFloat to start with */ 2774/* dfr is the decFloat to move toward */ 2775/* set is the context */ 2776/* returns result */ 2777/* */ 2778/* This is 754-1985 nextafter, as modified during revision (dropped */ 2779/* from 754-2008); status may be set unless the result is a normal */ 2780/* number. */ 2781/* ------------------------------------------------------------------ */ 2782decFloat * decFloatNextToward(decFloat *result, 2783 const decFloat *dfl, const decFloat *dfr, 2784 decContext *set) { 2785 decFloat delta; /* tiny increment or decrement */ 2786 decFloat pointone; /* 1e-1 */ 2787 uInt savestat; /* saves status */ 2788 enum rounding saveround; /* .. and mode */ 2789 uInt deltatop; /* top word for delta */ 2790 Int comp; /* work */ 2791 2792 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set); 2793 /* Both are numeric, so Invalid no longer a possibility */ 2794 comp=decNumCompare(dfl, dfr, 0); 2795 if (comp==0) return decFloatCopySign(result, dfl, dfr); /* equal */ 2796 /* unequal; do NextPlus or NextMinus but with different status rules */ 2797 2798 if (comp<0) { /* lhs<rhs, do NextPlus, see above for commentary */ 2799 if (DFISINF(dfl) && DFISSIGNED(dfl)) { /* -Infinity special case */ 2800 DFSETNMAX(result); 2801 DFWORD(result, 0)|=DECFLOAT_Sign; 2802 return result; 2803 } 2804 saveround=set->round; /* save mode */ 2805 set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */ 2806 deltatop=0; /* positive delta */ 2807 } 2808 else { /* lhs>rhs, do NextMinus, see above for commentary */ 2809 if (DFISINF(dfl) && !DFISSIGNED(dfl)) { /* +Infinity special case */ 2810 DFSETNMAX(result); 2811 return result; 2812 } 2813 saveround=set->round; /* save mode */ 2814 set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */ 2815 deltatop=DECFLOAT_Sign; /* negative delta */ 2816 } 2817 savestat=set->status; /* save status */ 2818 /* Here, Inexact is needed where appropriate (and hence Underflow, */ 2819 /* etc.). Therefore the tiny delta which is otherwise */ 2820 /* unrepresentable (see NextPlus and NextMinus) is constructed */ 2821 /* using the multiplication of FMA. */ 2822 decFloatZero(&delta); /* set up tiny delta */ 2823 DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */ 2824 DFWORD(&delta, 0)=deltatop; /* Sign + biased exponent=0 */ 2825 decFloatFromString(&pointone, "1E-1", set); /* set up multiplier */ 2826 decFloatFMA(result, &delta, &pointone, dfl, set); 2827 /* [Delta is truly tiny, so no need to correct sign of zero] */ 2828 /* use new status unless the result is normal */ 2829 if (decFloatIsNormal(result)) set->status=savestat; /* else goes forward */ 2830 set->round=saveround; /* restore mode */ 2831 return result; 2832 } /* decFloatNextToward */ 2833 2834/* ------------------------------------------------------------------ */ 2835/* decFloatOr -- logical digitwise OR of two decFloats */ 2836/* */ 2837/* result gets the result of ORing dfl and dfr */ 2838/* dfl is the first decFloat (lhs) */ 2839/* dfr is the second decFloat (rhs) */ 2840/* set is the context */ 2841/* returns result, which will be canonical with sign=0 */ 2842/* */ 2843/* The operands must be positive, finite with exponent q=0, and */ 2844/* comprise just zeros and ones; if not, Invalid operation results. */ 2845/* ------------------------------------------------------------------ */ 2846decFloat * decFloatOr(decFloat *result, 2847 const decFloat *dfl, const decFloat *dfr, 2848 decContext *set) { 2849 if (!DFISUINT01(dfl) || !DFISUINT01(dfr) 2850 || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set); 2851 /* the operands are positive finite integers (q=0) with just 0s and 1s */ 2852 #if DOUBLE 2853 DFWORD(result, 0)=ZEROWORD 2854 |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04009124); 2855 DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x49124491; 2856 #elif QUAD 2857 DFWORD(result, 0)=ZEROWORD 2858 |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04000912); 2859 DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x44912449; 2860 DFWORD(result, 2)=(DFWORD(dfl, 2) | DFWORD(dfr, 2))&0x12449124; 2861 DFWORD(result, 3)=(DFWORD(dfl, 3) | DFWORD(dfr, 3))&0x49124491; 2862 #endif 2863 return result; 2864 } /* decFloatOr */ 2865 2866/* ------------------------------------------------------------------ */ 2867/* decFloatPlus -- add value to 0, heeding NaNs, etc. */ 2868/* */ 2869/* result gets the canonicalized 0+df */ 2870/* df is the decFloat to plus */ 2871/* set is the context */ 2872/* returns result */ 2873/* */ 2874/* This has the same effect as 0+df where the exponent of the zero is */ 2875/* the same as that of df (if df is finite). */ 2876/* The effect is also the same as decFloatCopy except that NaNs */ 2877/* are handled normally (the sign of a NaN is not affected, and an */ 2878/* sNaN will signal), the result is canonical, and zero gets sign 0. */ 2879/* ------------------------------------------------------------------ */ 2880decFloat * decFloatPlus(decFloat *result, const decFloat *df, 2881 decContext *set) { 2882 if (DFISNAN(df)) return decNaNs(result, df, NULL, set); 2883 decCanonical(result, df); /* copy and check */ 2884 if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */ 2885 return result; 2886 } /* decFloatPlus */ 2887 2888/* ------------------------------------------------------------------ */ 2889/* decFloatQuantize -- quantize a decFloat */ 2890/* */ 2891/* result gets the result of quantizing dfl to match dfr */ 2892/* dfl is the first decFloat (lhs) */ 2893/* dfr is the second decFloat (rhs), which sets the exponent */ 2894/* set is the context */ 2895/* returns result */ 2896/* */ 2897/* Unless there is an error or the result is infinite, the exponent */ 2898/* of result is guaranteed to be the same as that of dfr. */ 2899/* ------------------------------------------------------------------ */ 2900decFloat * decFloatQuantize(decFloat *result, 2901 const decFloat *dfl, const decFloat *dfr, 2902 decContext *set) { 2903 Int explb, exprb; /* left and right biased exponents */ 2904 uByte *ulsd; /* local LSD pointer */ 2905 uByte *ub, *uc; /* work */ 2906 Int drop; /* .. */ 2907 uInt dpd; /* .. */ 2908 uInt encode; /* encoding accumulator */ 2909 uInt sourhil, sourhir; /* top words from source decFloats */ 2910 uInt uiwork; /* for macros */ 2911 #if QUAD 2912 uShort uswork; /* .. */ 2913 #endif 2914 /* the following buffer holds the coefficient for manipulation */ 2915 uByte buf[4+DECPMAX*3+2*QUAD]; /* + space for zeros to left or right */ 2916 #if DECTRACE 2917 bcdnum num; /* for trace displays */ 2918 #endif 2919 2920 /* Start decoding the arguments */ 2921 sourhil=DFWORD(dfl, 0); /* LHS top word */ 2922 explb=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */ 2923 sourhir=DFWORD(dfr, 0); /* RHS top word */ 2924 exprb=DECCOMBEXP[sourhir>>26]; 2925 2926 if (EXPISSPECIAL(explb | exprb)) { /* either is special? */ 2927 /* NaNs are handled as usual */ 2928 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set); 2929 /* one infinity but not both is bad */ 2930 if (DFISINF(dfl)!=DFISINF(dfr)) return decInvalid(result, set); 2931 /* both infinite; return canonical infinity with sign of LHS */ 2932 return decInfinity(result, dfl); 2933 } 2934 2935 /* Here when both arguments are finite */ 2936 /* complete extraction of the exponents [no need to unbias] */ 2937 explb+=GETECON(dfl); /* + continuation */ 2938 exprb+=GETECON(dfr); /* .. */ 2939 2940 /* calculate the number of digits to drop from the coefficient */ 2941 drop=exprb-explb; /* 0 if nothing to do */ 2942 if (drop==0) return decCanonical(result, dfl); /* return canonical */ 2943 2944 /* the coefficient is needed; lay it out into buf, offset so zeros */ 2945 /* can be added before or after as needed -- an extra heading is */ 2946 /* added so can safely pad Quad DECPMAX-1 zeros to the left by */ 2947 /* fours */ 2948 #define BUFOFF (buf+4+DECPMAX) 2949 GETCOEFF(dfl, BUFOFF); /* decode from decFloat */ 2950 /* [now the msd is at BUFOFF and the lsd is at BUFOFF+DECPMAX-1] */ 2951 2952 #if DECTRACE 2953 num.msd=BUFOFF; 2954 num.lsd=BUFOFF+DECPMAX-1; 2955 num.exponent=explb-DECBIAS; 2956 num.sign=sourhil & DECFLOAT_Sign; 2957 decShowNum(&num, "dfl"); 2958 #endif 2959 2960 if (drop>0) { /* [most common case] */ 2961 /* (this code is very similar to that in decFloatFinalize, but */ 2962 /* has many differences so is duplicated here -- so any changes */ 2963 /* may need to be made there, too) */ 2964 uByte *roundat; /* -> re-round digit */ 2965 uByte reround; /* reround value */ 2966 /* printf("Rounding; drop=%ld\n", (LI)drop); */ 2967 2968 /* there is at least one zero needed to the left, in all but one */ 2969 /* exceptional (all-nines) case, so place four zeros now; this is */ 2970 /* needed almost always and makes rounding all-nines by fours safe */ 2971 UBFROMUI(BUFOFF-4, 0); 2972 2973 /* Three cases here: */ 2974 /* 1. new LSD is in coefficient (almost always) */ 2975 /* 2. new LSD is digit to left of coefficient (so MSD is */ 2976 /* round-for-reround digit) */ 2977 /* 3. new LSD is to left of case 2 (whole coefficient is sticky) */ 2978 /* Note that leading zeros can safely be treated as useful digits */ 2979 2980 /* [duplicate check-stickies code to save a test] */ 2981 /* [by-digit check for stickies as runs of zeros are rare] */ 2982 if (drop<DECPMAX) { /* NB lengths not addresses */ 2983 roundat=BUFOFF+DECPMAX-drop; 2984 reround=*roundat; 2985 for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) { 2986 if (*ub!=0) { /* non-zero to be discarded */ 2987 reround=DECSTICKYTAB[reround]; /* apply sticky bit */ 2988 break; /* [remainder don't-care] */ 2989 } 2990 } /* check stickies */ 2991 ulsd=roundat-1; /* set LSD */ 2992 } 2993 else { /* edge case */ 2994 if (drop==DECPMAX) { 2995 roundat=BUFOFF; 2996 reround=*roundat; 2997 } 2998 else { 2999 roundat=BUFOFF-1; 3000 reround=0; 3001 } 3002 for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) { 3003 if (*ub!=0) { /* non-zero to be discarded */ 3004 reround=DECSTICKYTAB[reround]; /* apply sticky bit */ 3005 break; /* [remainder don't-care] */ 3006 } 3007 } /* check stickies */ 3008 *BUFOFF=0; /* make a coefficient of 0 */ 3009 ulsd=BUFOFF; /* .. at the MSD place */ 3010 } 3011 3012 if (reround!=0) { /* discarding non-zero */ 3013 uInt bump=0; 3014 set->status|=DEC_Inexact; 3015 3016 /* next decide whether to increment the coefficient */ 3017 if (set->round==DEC_ROUND_HALF_EVEN) { /* fastpath slowest case */ 3018 if (reround>5) bump=1; /* >0.5 goes up */ 3019 else if (reround==5) /* exactly 0.5000 .. */ 3020 bump=*ulsd & 0x01; /* .. up iff [new] lsd is odd */ 3021 } /* r-h-e */ 3022 else switch (set->round) { 3023 case DEC_ROUND_DOWN: { 3024 /* no change */ 3025 break;} /* r-d */ 3026 case DEC_ROUND_HALF_DOWN: { 3027 if (reround>5) bump=1; 3028 break;} /* r-h-d */ 3029 case DEC_ROUND_HALF_UP: { 3030 if (reround>=5) bump=1; 3031 break;} /* r-h-u */ 3032 case DEC_ROUND_UP: { 3033 if (reround>0) bump=1; 3034 break;} /* r-u */ 3035 case DEC_ROUND_CEILING: { 3036 /* same as _UP for positive numbers, and as _DOWN for negatives */ 3037 if (!(sourhil&DECFLOAT_Sign) && reround>0) bump=1; 3038 break;} /* r-c */ 3039 case DEC_ROUND_FLOOR: { 3040 /* same as _UP for negative numbers, and as _DOWN for positive */ 3041 /* [negative reround cannot occur on 0] */ 3042 if (sourhil&DECFLOAT_Sign && reround>0) bump=1; 3043 break;} /* r-f */ 3044 case DEC_ROUND_05UP: { 3045 if (reround>0) { /* anything out there is 'sticky' */ 3046 /* bump iff lsd=0 or 5; this cannot carry so it could be */ 3047 /* effected immediately with no bump -- but the code */ 3048 /* is clearer if this is done the same way as the others */ 3049 if (*ulsd==0 || *ulsd==5) bump=1; 3050 } 3051 break;} /* r-r */ 3052 default: { /* e.g., DEC_ROUND_MAX */ 3053 set->status|=DEC_Invalid_context; 3054 #if DECCHECK 3055 printf("Unknown rounding mode: %ld\n", (LI)set->round); 3056 #endif 3057 break;} 3058 } /* switch (not r-h-e) */ 3059 /* printf("ReRound: %ld bump: %ld\n", (LI)reround, (LI)bump); */ 3060 3061 if (bump!=0) { /* need increment */ 3062 /* increment the coefficient; this could give 1000... (after */ 3063 /* the all nines case) */ 3064 ub=ulsd; 3065 for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0); 3066 /* now at most 3 digits left to non-9 (usually just the one) */ 3067 for (; *ub==9; ub--) *ub=0; 3068 *ub+=1; 3069 /* [the all-nines case will have carried one digit to the */ 3070 /* left of the original MSD -- just where it is needed] */ 3071 } /* bump needed */ 3072 } /* inexact rounding */ 3073 3074 /* now clear zeros to the left so exactly DECPMAX digits will be */ 3075 /* available in the coefficent -- the first word to the left was */ 3076 /* cleared earlier for safe carry; now add any more needed */ 3077 if (drop>4) { 3078 UBFROMUI(BUFOFF-8, 0); /* must be at least 5 */ 3079 for (uc=BUFOFF-12; uc>ulsd-DECPMAX-3; uc-=4) UBFROMUI(uc, 0); 3080 } 3081 } /* need round (drop>0) */ 3082 3083 else { /* drop<0; padding with -drop digits is needed */ 3084 /* This is the case where an error can occur if the padded */ 3085 /* coefficient will not fit; checking for this can be done in the */ 3086 /* same loop as padding for zeros if the no-hope and zero cases */ 3087 /* are checked first */ 3088 if (-drop>DECPMAX-1) { /* cannot fit unless 0 */ 3089 if (!ISCOEFFZERO(BUFOFF)) return decInvalid(result, set); 3090 /* a zero can have any exponent; just drop through and use it */ 3091 ulsd=BUFOFF+DECPMAX-1; 3092 } 3093 else { /* padding will fit (but may still be too long) */ 3094 /* final-word mask depends on endianess */ 3095 #if DECLITEND 3096 static const uInt dmask[]={0, 0x000000ff, 0x0000ffff, 0x00ffffff}; 3097 #else 3098 static const uInt dmask[]={0, 0xff000000, 0xffff0000, 0xffffff00}; 3099 #endif 3100 /* note that here zeros to the right are added by fours, so in */ 3101 /* the Quad case this could write 36 zeros if the coefficient has */ 3102 /* fewer than three significant digits (hence the +2*QUAD for buf) */ 3103 for (uc=BUFOFF+DECPMAX;; uc+=4) { 3104 UBFROMUI(uc, 0); 3105 if (UBTOUI(uc-DECPMAX)!=0) { /* could be bad */ 3106 /* if all four digits should be zero, definitely bad */ 3107 if (uc<=BUFOFF+DECPMAX+(-drop)-4) 3108 return decInvalid(result, set); 3109 /* must be a 1- to 3-digit sequence; check more carefully */ 3110 if ((UBTOUI(uc-DECPMAX)&dmask[(-drop)%4])!=0) 3111 return decInvalid(result, set); 3112 break; /* no need for loop end test */ 3113 } 3114 if (uc>=BUFOFF+DECPMAX+(-drop)-4) break; /* done */ 3115 } 3116 ulsd=BUFOFF+DECPMAX+(-drop)-1; 3117 } /* pad and check leading zeros */ 3118 } /* drop<0 */ 3119 3120 #if DECTRACE 3121 num.msd=ulsd-DECPMAX+1; 3122 num.lsd=ulsd; 3123 num.exponent=explb-DECBIAS; 3124 num.sign=sourhil & DECFLOAT_Sign; 3125 decShowNum(&num, "res"); 3126 #endif 3127 3128 /*------------------------------------------------------------------*/ 3129 /* At this point the result is DECPMAX digits, ending at ulsd, so */ 3130 /* fits the encoding exactly; there is no possibility of error */ 3131 /*------------------------------------------------------------------*/ 3132 encode=((exprb>>DECECONL)<<4) + *(ulsd-DECPMAX+1); /* make index */ 3133 encode=DECCOMBFROM[encode]; /* indexed by (0-2)*16+msd */ 3134 /* the exponent continuation can be extracted from the original RHS */ 3135 encode|=sourhir & ECONMASK; 3136 encode|=sourhil&DECFLOAT_Sign; /* add the sign from LHS */ 3137 3138 /* finally encode the coefficient */ 3139 /* private macro to encode a declet; this version can be used */ 3140 /* because all coefficient digits exist */ 3141 #define getDPD3q(dpd, n) ub=ulsd-(3*(n))-2; \ 3142 dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)]; 3143 3144 #if DOUBLE 3145 getDPD3q(dpd, 4); encode|=dpd<<8; 3146 getDPD3q(dpd, 3); encode|=dpd>>2; 3147 DFWORD(result, 0)=encode; 3148 encode=dpd<<30; 3149 getDPD3q(dpd, 2); encode|=dpd<<20; 3150 getDPD3q(dpd, 1); encode|=dpd<<10; 3151 getDPD3q(dpd, 0); encode|=dpd; 3152 DFWORD(result, 1)=encode; 3153 3154 #elif QUAD 3155 getDPD3q(dpd,10); encode|=dpd<<4; 3156 getDPD3q(dpd, 9); encode|=dpd>>6; 3157 DFWORD(result, 0)=encode; 3158 encode=dpd<<26; 3159 getDPD3q(dpd, 8); encode|=dpd<<16; 3160 getDPD3q(dpd, 7); encode|=dpd<<6; 3161 getDPD3q(dpd, 6); encode|=dpd>>4; 3162 DFWORD(result, 1)=encode; 3163 encode=dpd<<28; 3164 getDPD3q(dpd, 5); encode|=dpd<<18; 3165 getDPD3q(dpd, 4); encode|=dpd<<8; 3166 getDPD3q(dpd, 3); encode|=dpd>>2; 3167 DFWORD(result, 2)=encode; 3168 encode=dpd<<30; 3169 getDPD3q(dpd, 2); encode|=dpd<<20; 3170 getDPD3q(dpd, 1); encode|=dpd<<10; 3171 getDPD3q(dpd, 0); encode|=dpd; 3172 DFWORD(result, 3)=encode; 3173 #endif 3174 return result; 3175 } /* decFloatQuantize */ 3176 3177/* ------------------------------------------------------------------ */ 3178/* decFloatReduce -- reduce finite coefficient to minimum length */ 3179/* */ 3180/* result gets the reduced decFloat */ 3181/* df is the source decFloat */ 3182/* set is the context */ 3183/* returns result, which will be canonical */ 3184/* */ 3185/* This removes all possible trailing zeros from the coefficient; */ 3186/* some may remain when the number is very close to Nmax. */ 3187/* Special values are unchanged and no status is set unless df=sNaN. */ 3188/* Reduced zero has an exponent q=0. */ 3189/* ------------------------------------------------------------------ */ 3190decFloat * decFloatReduce(decFloat *result, const decFloat *df, 3191 decContext *set) { 3192 bcdnum num; /* work */ 3193 uByte buf[DECPMAX], *ub; /* coefficient and pointer */ 3194 if (df!=result) *result=*df; /* copy, if needed */ 3195 if (DFISNAN(df)) return decNaNs(result, df, NULL, set); /* sNaN */ 3196 /* zeros and infinites propagate too */ 3197 if (DFISINF(df)) return decInfinity(result, df); /* canonical */ 3198 if (DFISZERO(df)) { 3199 uInt sign=DFWORD(df, 0)&DECFLOAT_Sign; 3200 decFloatZero(result); 3201 DFWORD(result, 0)|=sign; 3202 return result; /* exponent dropped, sign OK */ 3203 } 3204 /* non-zero finite */ 3205 GETCOEFF(df, buf); 3206 ub=buf+DECPMAX-1; /* -> lsd */ 3207 if (*ub) return result; /* no trailing zeros */ 3208 for (ub--; *ub==0;) ub--; /* terminates because non-zero */ 3209 /* *ub is the first non-zero from the right */ 3210 num.sign=DFWORD(df, 0)&DECFLOAT_Sign; /* set up number... */ 3211 num.exponent=GETEXPUN(df)+(Int)(buf+DECPMAX-1-ub); /* adjusted exponent */ 3212 num.msd=buf; 3213 num.lsd=ub; 3214 return decFinalize(result, &num, set); 3215 } /* decFloatReduce */ 3216 3217/* ------------------------------------------------------------------ */ 3218/* decFloatRemainder -- integer divide and return remainder */ 3219/* */ 3220/* result gets the remainder of dividing dfl by dfr: */ 3221/* dfl is the first decFloat (lhs) */ 3222/* dfr is the second decFloat (rhs) */ 3223/* set is the context */ 3224/* returns result */ 3225/* */ 3226/* ------------------------------------------------------------------ */ 3227decFloat * decFloatRemainder(decFloat *result, 3228 const decFloat *dfl, const decFloat *dfr, 3229 decContext *set) { 3230 return decDivide(result, dfl, dfr, set, REMAINDER); 3231 } /* decFloatRemainder */ 3232 3233/* ------------------------------------------------------------------ */ 3234/* decFloatRemainderNear -- integer divide to nearest and remainder */ 3235/* */ 3236/* result gets the remainder of dividing dfl by dfr: */ 3237/* dfl is the first decFloat (lhs) */ 3238/* dfr is the second decFloat (rhs) */ 3239/* set is the context */ 3240/* returns result */ 3241/* */ 3242/* This is the IEEE remainder, where the nearest integer is used. */ 3243/* ------------------------------------------------------------------ */ 3244decFloat * decFloatRemainderNear(decFloat *result, 3245 const decFloat *dfl, const decFloat *dfr, 3246 decContext *set) { 3247 return decDivide(result, dfl, dfr, set, REMNEAR); 3248 } /* decFloatRemainderNear */ 3249 3250/* ------------------------------------------------------------------ */ 3251/* decFloatRotate -- rotate the coefficient of a decFloat left/right */ 3252/* */ 3253/* result gets the result of rotating dfl */ 3254/* dfl is the source decFloat to rotate */ 3255/* dfr is the count of digits to rotate, an integer (with q=0) */ 3256/* set is the context */ 3257/* returns result */ 3258/* */ 3259/* The digits of the coefficient of dfl are rotated to the left (if */ 3260/* dfr is positive) or to the right (if dfr is negative) without */ 3261/* adjusting the exponent or the sign of dfl. */ 3262/* */ 3263/* dfr must be in the range -DECPMAX through +DECPMAX. */ 3264/* NaNs are propagated as usual. An infinite dfl is unaffected (but */ 3265/* dfr must be valid). No status is set unless dfr is invalid or an */ 3266/* operand is an sNaN. The result is canonical. */ 3267/* ------------------------------------------------------------------ */ 3268#define PHALF (ROUNDUP(DECPMAX/2, 4)) /* half length, rounded up */ 3269decFloat * decFloatRotate(decFloat *result, 3270 const decFloat *dfl, const decFloat *dfr, 3271 decContext *set) { 3272 Int rotate; /* dfr as an Int */ 3273 uByte buf[DECPMAX+PHALF]; /* coefficient + half */ 3274 uInt digits, savestat; /* work */ 3275 bcdnum num; /* .. */ 3276 uByte *ub; /* .. */ 3277 3278 if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set); 3279 if (!DFISINT(dfr)) return decInvalid(result, set); 3280 digits=decFloatDigits(dfr); /* calculate digits */ 3281 if (digits>2) return decInvalid(result, set); /* definitely out of range */ 3282 rotate=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */ 3283 if (rotate>DECPMAX) return decInvalid(result, set); /* too big */ 3284 /* [from here on no error or status change is possible] */ 3285 if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */ 3286 /* handle no-rotate cases */ 3287 if (rotate==0 || rotate==DECPMAX) return decCanonical(result, dfl); 3288 /* a real rotate is needed: 0 < rotate < DECPMAX */ 3289 /* reduce the rotation to no more than half to reduce copying later */ 3290 /* (for QUAD in fact half + 2 digits) */ 3291 if (DFISSIGNED(dfr)) rotate=-rotate; 3292 if (abs(rotate)>PHALF) { 3293 if (rotate<0) rotate=DECPMAX+rotate; 3294 else rotate=rotate-DECPMAX; 3295 } 3296 /* now lay out the coefficient, leaving room to the right or the */ 3297 /* left depending on the direction of rotation */ 3298 ub=buf; 3299 if (rotate<0) ub+=PHALF; /* rotate right, so space to left */ 3300 GETCOEFF(dfl, ub); 3301 /* copy half the digits to left or right, and set num.msd */ 3302 if (rotate<0) { 3303 memcpy(buf, buf+DECPMAX, PHALF); 3304 num.msd=buf+PHALF+rotate; 3305 } 3306 else { 3307 memcpy(buf+DECPMAX, buf, PHALF); 3308 num.msd=buf+rotate; 3309 } 3310 /* fill in rest of num */ 3311 num.lsd=num.msd+DECPMAX-1; 3312 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign; 3313 num.exponent=GETEXPUN(dfl); 3314 savestat=set->status; /* record */ 3315 decFinalize(result, &num, set); 3316 set->status=savestat; /* restore */ 3317 return result; 3318 } /* decFloatRotate */ 3319 3320/* ------------------------------------------------------------------ */ 3321/* decFloatSameQuantum -- test decFloats for same quantum */ 3322/* */ 3323/* dfl is the first decFloat (lhs) */ 3324/* dfr is the second decFloat (rhs) */ 3325/* returns 1 if the operands have the same quantum, 0 otherwise */ 3326/* */ 3327/* No error is possible and no status results. */ 3328/* ------------------------------------------------------------------ */ 3329uInt decFloatSameQuantum(const decFloat *dfl, const decFloat *dfr) { 3330 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { 3331 if (DFISNAN(dfl) && DFISNAN(dfr)) return 1; 3332 if (DFISINF(dfl) && DFISINF(dfr)) return 1; 3333 return 0; /* any other special mixture gives false */ 3334 } 3335 if (GETEXP(dfl)==GETEXP(dfr)) return 1; /* biased exponents match */ 3336 return 0; 3337 } /* decFloatSameQuantum */ 3338 3339/* ------------------------------------------------------------------ */ 3340/* decFloatScaleB -- multiply by a power of 10, as per 754 */ 3341/* */ 3342/* result gets the result of the operation */ 3343/* dfl is the first decFloat (lhs) */ 3344/* dfr is the second decFloat (rhs), am integer (with q=0) */ 3345/* set is the context */ 3346/* returns result */ 3347/* */ 3348/* This computes result=dfl x 10**dfr where dfr is an integer in the */ 3349/* range +/-2*(emax+pmax), typically resulting from LogB. */ 3350/* Underflow and Overflow (with Inexact) may occur. NaNs propagate */ 3351/* as usual. */ 3352/* ------------------------------------------------------------------ */ 3353#define SCALEBMAX 2*(DECEMAX+DECPMAX) /* D=800, Q=12356 */ 3354decFloat * decFloatScaleB(decFloat *result, 3355 const decFloat *dfl, const decFloat *dfr, 3356 decContext *set) { 3357 uInt digits; /* work */ 3358 Int expr; /* dfr as an Int */ 3359 3360 if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set); 3361 if (!DFISINT(dfr)) return decInvalid(result, set); 3362 digits=decFloatDigits(dfr); /* calculate digits */ 3363 3364 #if DOUBLE 3365 if (digits>3) return decInvalid(result, set); /* definitely out of range */ 3366 expr=DPD2BIN[DFWORD(dfr, 1)&0x3ff]; /* must be in bottom declet */ 3367 #elif QUAD 3368 if (digits>5) return decInvalid(result, set); /* definitely out of range */ 3369 expr=DPD2BIN[DFWORD(dfr, 3)&0x3ff] /* in bottom 2 declets .. */ 3370 +DPD2BIN[(DFWORD(dfr, 3)>>10)&0x3ff]*1000; /* .. */ 3371 #endif 3372 if (expr>SCALEBMAX) return decInvalid(result, set); /* oops */ 3373 /* [from now on no error possible] */ 3374 if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */ 3375 if (DFISSIGNED(dfr)) expr=-expr; 3376 /* dfl is finite and expr is valid */ 3377 *result=*dfl; /* copy to target */ 3378 return decFloatSetExponent(result, set, GETEXPUN(result)+expr); 3379 } /* decFloatScaleB */ 3380 3381/* ------------------------------------------------------------------ */ 3382/* decFloatShift -- shift the coefficient of a decFloat left or right */ 3383/* */ 3384/* result gets the result of shifting dfl */ 3385/* dfl is the source decFloat to shift */ 3386/* dfr is the count of digits to shift, an integer (with q=0) */ 3387/* set is the context */ 3388/* returns result */ 3389/* */ 3390/* The digits of the coefficient of dfl are shifted to the left (if */ 3391/* dfr is positive) or to the right (if dfr is negative) without */ 3392/* adjusting the exponent or the sign of dfl. */ 3393/* */ 3394/* dfr must be in the range -DECPMAX through +DECPMAX. */ 3395/* NaNs are propagated as usual. An infinite dfl is unaffected (but */ 3396/* dfr must be valid). No status is set unless dfr is invalid or an */ 3397/* operand is an sNaN. The result is canonical. */ 3398/* ------------------------------------------------------------------ */ 3399decFloat * decFloatShift(decFloat *result, 3400 const decFloat *dfl, const decFloat *dfr, 3401 decContext *set) { 3402 Int shift; /* dfr as an Int */ 3403 uByte buf[DECPMAX*2]; /* coefficient + padding */ 3404 uInt digits, savestat; /* work */ 3405 bcdnum num; /* .. */ 3406 uInt uiwork; /* for macros */ 3407 3408 if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set); 3409 if (!DFISINT(dfr)) return decInvalid(result, set); 3410 digits=decFloatDigits(dfr); /* calculate digits */ 3411 if (digits>2) return decInvalid(result, set); /* definitely out of range */ 3412 shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */ 3413 if (shift>DECPMAX) return decInvalid(result, set); /* too big */ 3414 /* [from here on no error or status change is possible] */ 3415 3416 if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */ 3417 /* handle no-shift and all-shift (clear to zero) cases */ 3418 if (shift==0) return decCanonical(result, dfl); 3419 if (shift==DECPMAX) { /* zero with sign */ 3420 uByte sign=(uByte)(DFBYTE(dfl, 0)&0x80); /* save sign bit */ 3421 decFloatZero(result); /* make +0 */ 3422 DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* and set sign */ 3423 /* [cannot safely use CopySign] */ 3424 return result; 3425 } 3426 /* a real shift is needed: 0 < shift < DECPMAX */ 3427 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign; 3428 num.exponent=GETEXPUN(dfl); 3429 num.msd=buf; 3430 GETCOEFF(dfl, buf); 3431 if (DFISSIGNED(dfr)) { /* shift right */ 3432 /* edge cases are taken care of, so this is easy */ 3433 num.lsd=buf+DECPMAX-shift-1; 3434 } 3435 else { /* shift left -- zero padding needed to right */ 3436 UBFROMUI(buf+DECPMAX, 0); /* 8 will handle most cases */ 3437 UBFROMUI(buf+DECPMAX+4, 0); /* .. */ 3438 if (shift>8) memset(buf+DECPMAX+8, 0, 8+QUAD*18); /* all other cases */ 3439 num.msd+=shift; 3440 num.lsd=num.msd+DECPMAX-1; 3441 } 3442 savestat=set->status; /* record */ 3443 decFinalize(result, &num, set); 3444 set->status=savestat; /* restore */ 3445 return result; 3446 } /* decFloatShift */ 3447 3448/* ------------------------------------------------------------------ */ 3449/* decFloatSubtract -- subtract a decFloat from another */ 3450/* */ 3451/* result gets the result of subtracting dfr from dfl: */ 3452/* dfl is the first decFloat (lhs) */ 3453/* dfr is the second decFloat (rhs) */ 3454/* set is the context */ 3455/* returns result */ 3456/* */ 3457/* ------------------------------------------------------------------ */ 3458decFloat * decFloatSubtract(decFloat *result, 3459 const decFloat *dfl, const decFloat *dfr, 3460 decContext *set) { 3461 decFloat temp; 3462 /* NaNs must propagate without sign change */ 3463 if (DFISNAN(dfr)) return decFloatAdd(result, dfl, dfr, set); 3464 temp=*dfr; /* make a copy */ 3465 DFBYTE(&temp, 0)^=0x80; /* flip sign */ 3466 return decFloatAdd(result, dfl, &temp, set); /* and add to the lhs */ 3467 } /* decFloatSubtract */ 3468 3469/* ------------------------------------------------------------------ */ 3470/* decFloatToInt -- round to 32-bit binary integer (4 flavours) */ 3471/* */ 3472/* df is the decFloat to round */ 3473/* set is the context */ 3474/* round is the rounding mode to use */ 3475/* returns a uInt or an Int, rounded according to the name */ 3476/* */ 3477/* Invalid will always be signaled if df is a NaN, is Infinite, or is */ 3478/* outside the range of the target; Inexact will not be signaled for */ 3479/* simple rounding unless 'Exact' appears in the name. */ 3480/* ------------------------------------------------------------------ */ 3481uInt decFloatToUInt32(const decFloat *df, decContext *set, 3482 enum rounding round) { 3483 return decToInt32(df, set, round, 0, 1);} 3484 3485uInt decFloatToUInt32Exact(const decFloat *df, decContext *set, 3486 enum rounding round) { 3487 return decToInt32(df, set, round, 1, 1);} 3488 3489Int decFloatToInt32(const decFloat *df, decContext *set, 3490 enum rounding round) { 3491 return (Int)decToInt32(df, set, round, 0, 0);} 3492 3493Int decFloatToInt32Exact(const decFloat *df, decContext *set, 3494 enum rounding round) { 3495 return (Int)decToInt32(df, set, round, 1, 0);} 3496 3497/* ------------------------------------------------------------------ */ 3498/* decFloatToIntegral -- round to integral value (two flavours) */ 3499/* */ 3500/* result gets the result */ 3501/* df is the decFloat to round */ 3502/* set is the context */ 3503/* round is the rounding mode to use */ 3504/* returns result */ 3505/* */ 3506/* No exceptions, even Inexact, are raised except for sNaN input, or */ 3507/* if 'Exact' appears in the name. */ 3508/* ------------------------------------------------------------------ */ 3509decFloat * decFloatToIntegralValue(decFloat *result, const decFloat *df, 3510 decContext *set, enum rounding round) { 3511 return decToIntegral(result, df, set, round, 0);} 3512 3513decFloat * decFloatToIntegralExact(decFloat *result, const decFloat *df, 3514 decContext *set) { 3515 return decToIntegral(result, df, set, set->round, 1);} 3516 3517/* ------------------------------------------------------------------ */ 3518/* decFloatXor -- logical digitwise XOR of two decFloats */ 3519/* */ 3520/* result gets the result of XORing dfl and dfr */ 3521/* dfl is the first decFloat (lhs) */ 3522/* dfr is the second decFloat (rhs) */ 3523/* set is the context */ 3524/* returns result, which will be canonical with sign=0 */ 3525/* */ 3526/* The operands must be positive, finite with exponent q=0, and */ 3527/* comprise just zeros and ones; if not, Invalid operation results. */ 3528/* ------------------------------------------------------------------ */ 3529decFloat * decFloatXor(decFloat *result, 3530 const decFloat *dfl, const decFloat *dfr, 3531 decContext *set) { 3532 if (!DFISUINT01(dfl) || !DFISUINT01(dfr) 3533 || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set); 3534 /* the operands are positive finite integers (q=0) with just 0s and 1s */ 3535 #if DOUBLE 3536 DFWORD(result, 0)=ZEROWORD 3537 |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04009124); 3538 DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x49124491; 3539 #elif QUAD 3540 DFWORD(result, 0)=ZEROWORD 3541 |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04000912); 3542 DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x44912449; 3543 DFWORD(result, 2)=(DFWORD(dfl, 2) ^ DFWORD(dfr, 2))&0x12449124; 3544 DFWORD(result, 3)=(DFWORD(dfl, 3) ^ DFWORD(dfr, 3))&0x49124491; 3545 #endif 3546 return result; 3547 } /* decFloatXor */ 3548 3549/* ------------------------------------------------------------------ */ 3550/* decInvalid -- set Invalid_operation result */ 3551/* */ 3552/* result gets a canonical NaN */ 3553/* set is the context */ 3554/* returns result */ 3555/* */ 3556/* status has Invalid_operation added */ 3557/* ------------------------------------------------------------------ */ 3558static decFloat *decInvalid(decFloat *result, decContext *set) { 3559 decFloatZero(result); 3560 DFWORD(result, 0)=DECFLOAT_qNaN; 3561 set->status|=DEC_Invalid_operation; 3562 return result; 3563 } /* decInvalid */ 3564 3565/* ------------------------------------------------------------------ */ 3566/* decInfinity -- set canonical Infinity with sign from a decFloat */ 3567/* */ 3568/* result gets a canonical Infinity */ 3569/* df is source decFloat (only the sign is used) */ 3570/* returns result */ 3571/* */ 3572/* df may be the same as result */ 3573/* ------------------------------------------------------------------ */ 3574static decFloat *decInfinity(decFloat *result, const decFloat *df) { 3575 uInt sign=DFWORD(df, 0); /* save source signword */ 3576 decFloatZero(result); /* clear everything */ 3577 DFWORD(result, 0)=DECFLOAT_Inf | (sign & DECFLOAT_Sign); 3578 return result; 3579 } /* decInfinity */ 3580 3581/* ------------------------------------------------------------------ */ 3582/* decNaNs -- handle NaN argument(s) */ 3583/* */ 3584/* result gets the result of handling dfl and dfr, one or both of */ 3585/* which is a NaN */ 3586/* dfl is the first decFloat (lhs) */ 3587/* dfr is the second decFloat (rhs) -- may be NULL for a single- */ 3588/* operand operation */ 3589/* set is the context */ 3590/* returns result */ 3591/* */ 3592/* Called when one or both operands is a NaN, and propagates the */ 3593/* appropriate result to res. When an sNaN is found, it is changed */ 3594/* to a qNaN and Invalid operation is set. */ 3595/* ------------------------------------------------------------------ */ 3596static decFloat *decNaNs(decFloat *result, 3597 const decFloat *dfl, const decFloat *dfr, 3598 decContext *set) { 3599 /* handle sNaNs first */ 3600 if (dfr!=NULL && DFISSNAN(dfr) && !DFISSNAN(dfl)) dfl=dfr; /* use RHS */ 3601 if (DFISSNAN(dfl)) { 3602 decCanonical(result, dfl); /* propagate canonical sNaN */ 3603 DFWORD(result, 0)&=~(DECFLOAT_qNaN ^ DECFLOAT_sNaN); /* quiet */ 3604 set->status|=DEC_Invalid_operation; 3605 return result; 3606 } 3607 /* one or both is a quiet NaN */ 3608 if (!DFISNAN(dfl)) dfl=dfr; /* RHS must be NaN, use it */ 3609 return decCanonical(result, dfl); /* propagate canonical qNaN */ 3610 } /* decNaNs */ 3611 3612/* ------------------------------------------------------------------ */ 3613/* decNumCompare -- numeric comparison of two decFloats */ 3614/* */ 3615/* dfl is the left-hand decFloat, which is not a NaN */ 3616/* dfr is the right-hand decFloat, which is not a NaN */ 3617/* tot is 1 for total order compare, 0 for simple numeric */ 3618/* returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr */ 3619/* */ 3620/* No error is possible; status and mode are unchanged. */ 3621/* ------------------------------------------------------------------ */ 3622static Int decNumCompare(const decFloat *dfl, const decFloat *dfr, Flag tot) { 3623 Int sigl, sigr; /* LHS and RHS non-0 signums */ 3624 Int shift; /* shift needed to align operands */ 3625 uByte *ub, *uc; /* work */ 3626 uInt uiwork; /* for macros */ 3627 /* buffers +2 if Quad (36 digits), need double plus 4 for safe padding */ 3628 uByte bufl[DECPMAX*2+QUAD*2+4]; /* for LHS coefficient + padding */ 3629 uByte bufr[DECPMAX*2+QUAD*2+4]; /* for RHS coefficient + padding */ 3630 3631 sigl=1; 3632 if (DFISSIGNED(dfl)) { 3633 if (!DFISSIGNED(dfr)) { /* -LHS +RHS */ 3634 if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0; 3635 return -1; /* RHS wins */ 3636 } 3637 sigl=-1; 3638 } 3639 if (DFISSIGNED(dfr)) { 3640 if (!DFISSIGNED(dfl)) { /* +LHS -RHS */ 3641 if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0; 3642 return +1; /* LHS wins */ 3643 } 3644 } 3645 3646 /* signs are the same; operand(s) could be zero */ 3647 sigr=-sigl; /* sign to return if abs(RHS) wins */ 3648 3649 if (DFISINF(dfl)) { 3650 if (DFISINF(dfr)) return 0; /* both infinite & same sign */ 3651 return sigl; /* inf > n */ 3652 } 3653 if (DFISINF(dfr)) return sigr; /* n < inf [dfl is finite] */ 3654 3655 /* here, both are same sign and finite; calculate their offset */ 3656 shift=GETEXP(dfl)-GETEXP(dfr); /* [0 means aligned] */ 3657 /* [bias can be ignored -- the absolute exponent is not relevant] */ 3658 3659 if (DFISZERO(dfl)) { 3660 if (!DFISZERO(dfr)) return sigr; /* LHS=0, RHS!=0 */ 3661 /* both are zero, return 0 if both same exponent or numeric compare */ 3662 if (shift==0 || !tot) return 0; 3663 if (shift>0) return sigl; 3664 return sigr; /* [shift<0] */ 3665 } 3666 else { /* LHS!=0 */ 3667 if (DFISZERO(dfr)) return sigl; /* LHS!=0, RHS=0 */ 3668 } 3669 /* both are known to be non-zero at this point */ 3670 3671 /* if the exponents are so different that the coefficients do not */ 3672 /* overlap (by even one digit) then a full comparison is not needed */ 3673 if (abs(shift)>=DECPMAX) { /* no overlap */ 3674 /* coefficients are known to be non-zero */ 3675 if (shift>0) return sigl; 3676 return sigr; /* [shift<0] */ 3677 } 3678 3679 /* decode the coefficients */ 3680 /* (shift both right two if Quad to make a multiple of four) */ 3681 #if QUAD 3682 UBFROMUI(bufl, 0); 3683 UBFROMUI(bufr, 0); 3684 #endif 3685 GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */ 3686 GETCOEFF(dfr, bufr+QUAD*2); /* .. */ 3687 if (shift==0) { /* aligned; common and easy */ 3688 /* all multiples of four, here */ 3689 for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) { 3690 uInt ui=UBTOUI(ub); 3691 if (ui==UBTOUI(uc)) continue; /* so far so same */ 3692 /* about to find a winner; go by bytes in case little-endian */ 3693 for (;; ub++, uc++) { 3694 if (*ub>*uc) return sigl; /* difference found */ 3695 if (*ub<*uc) return sigr; /* .. */ 3696 } 3697 } 3698 } /* aligned */ 3699 else if (shift>0) { /* lhs to left */ 3700 ub=bufl; /* RHS pointer */ 3701 /* pad bufl so right-aligned; most shifts will fit in 8 */ 3702 UBFROMUI(bufl+DECPMAX+QUAD*2, 0); /* add eight zeros */ 3703 UBFROMUI(bufl+DECPMAX+QUAD*2+4, 0); /* .. */ 3704 if (shift>8) { 3705 /* more than eight; fill the rest, and also worth doing the */ 3706 /* lead-in by fours */ 3707 uByte *up; /* work */ 3708 uByte *upend=bufl+DECPMAX+QUAD*2+shift; 3709 for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0); 3710 /* [pads up to 36 in all for Quad] */ 3711 for (;; ub+=4) { 3712 if (UBTOUI(ub)!=0) return sigl; 3713 if (ub+4>bufl+shift-4) break; 3714 } 3715 } 3716 /* check remaining leading digits */ 3717 for (; ub<bufl+shift; ub++) if (*ub!=0) return sigl; 3718 /* now start the overlapped part; bufl has been padded, so the */ 3719 /* comparison can go for the full length of bufr, which is a */ 3720 /* multiple of 4 bytes */ 3721 for (uc=bufr; ; uc+=4, ub+=4) { 3722 uInt ui=UBTOUI(ub); 3723 if (ui!=UBTOUI(uc)) { /* mismatch found */ 3724 for (;; uc++, ub++) { /* check from left [little-endian?] */ 3725 if (*ub>*uc) return sigl; /* difference found */ 3726 if (*ub<*uc) return sigr; /* .. */ 3727 } 3728 } /* mismatch */ 3729 if (uc==bufr+QUAD*2+DECPMAX-4) break; /* all checked */ 3730 } 3731 } /* shift>0 */ 3732 3733 else { /* shift<0) .. RHS is to left of LHS; mirror shift>0 */ 3734 uc=bufr; /* RHS pointer */ 3735 /* pad bufr so right-aligned; most shifts will fit in 8 */ 3736 UBFROMUI(bufr+DECPMAX+QUAD*2, 0); /* add eight zeros */ 3737 UBFROMUI(bufr+DECPMAX+QUAD*2+4, 0); /* .. */ 3738 if (shift<-8) { 3739 /* more than eight; fill the rest, and also worth doing the */ 3740 /* lead-in by fours */ 3741 uByte *up; /* work */ 3742 uByte *upend=bufr+DECPMAX+QUAD*2-shift; 3743 for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0); 3744 /* [pads up to 36 in all for Quad] */ 3745 for (;; uc+=4) { 3746 if (UBTOUI(uc)!=0) return sigr; 3747 if (uc+4>bufr-shift-4) break; 3748 } 3749 } 3750 /* check remaining leading digits */ 3751 for (; uc<bufr-shift; uc++) if (*uc!=0) return sigr; 3752 /* now start the overlapped part; bufr has been padded, so the */ 3753 /* comparison can go for the full length of bufl, which is a */ 3754 /* multiple of 4 bytes */ 3755 for (ub=bufl; ; ub+=4, uc+=4) { 3756 uInt ui=UBTOUI(ub); 3757 if (ui!=UBTOUI(uc)) { /* mismatch found */ 3758 for (;; ub++, uc++) { /* check from left [little-endian?] */ 3759 if (*ub>*uc) return sigl; /* difference found */ 3760 if (*ub<*uc) return sigr; /* .. */ 3761 } 3762 } /* mismatch */ 3763 if (ub==bufl+QUAD*2+DECPMAX-4) break; /* all checked */ 3764 } 3765 } /* shift<0 */ 3766 3767 /* Here when compare equal */ 3768 if (!tot) return 0; /* numerically equal */ 3769 /* total ordering .. exponent matters */ 3770 if (shift>0) return sigl; /* total order by exponent */ 3771 if (shift<0) return sigr; /* .. */ 3772 return 0; 3773 } /* decNumCompare */ 3774 3775/* ------------------------------------------------------------------ */ 3776/* decToInt32 -- local routine to effect ToInteger conversions */ 3777/* */ 3778/* df is the decFloat to convert */ 3779/* set is the context */ 3780/* rmode is the rounding mode to use */ 3781/* exact is 1 if Inexact should be signalled */ 3782/* unsign is 1 if the result a uInt, 0 if an Int (cast to uInt) */ 3783/* returns 32-bit result as a uInt */ 3784/* */ 3785/* Invalid is set is df is a NaN, is infinite, or is out-of-range; in */ 3786/* these cases 0 is returned. */ 3787/* ------------------------------------------------------------------ */ 3788static uInt decToInt32(const decFloat *df, decContext *set, 3789 enum rounding rmode, Flag exact, Flag unsign) { 3790 Int exp; /* exponent */ 3791 uInt sourhi, sourpen, sourlo; /* top word from source decFloat .. */ 3792 uInt hi, lo; /* .. penultimate, least, etc. */ 3793 decFloat zero, result; /* work */ 3794 Int i; /* .. */ 3795 3796 /* Start decoding the argument */ 3797 sourhi=DFWORD(df, 0); /* top word */ 3798 exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */ 3799 if (EXPISSPECIAL(exp)) { /* is special? */ 3800 set->status|=DEC_Invalid_operation; /* signal */ 3801 return 0; 3802 } 3803 3804 /* Here when the argument is finite */ 3805 if (GETEXPUN(df)==0) result=*df; /* already a true integer */ 3806 else { /* need to round to integer */ 3807 enum rounding saveround; /* saver */ 3808 uInt savestatus; /* .. */ 3809 saveround=set->round; /* save rounding mode .. */ 3810 savestatus=set->status; /* .. and status */ 3811 set->round=rmode; /* set mode */ 3812 decFloatZero(&zero); /* make 0E+0 */ 3813 set->status=0; /* clear */ 3814 decFloatQuantize(&result, df, &zero, set); /* [this may fail] */ 3815 set->round=saveround; /* restore rounding mode .. */ 3816 if (exact) set->status|=savestatus; /* include Inexact */ 3817 else set->status=savestatus; /* .. or just original status */ 3818 } 3819 3820 /* only the last four declets of the coefficient can contain */ 3821 /* non-zero; check for others (and also NaN or Infinity from the */ 3822 /* Quantize) first (see DFISZERO for explanation): */ 3823 /* decFloatShow(&result, "sofar"); */ 3824 #if DOUBLE 3825 if ((DFWORD(&result, 0)&0x1c03ff00)!=0 3826 || (DFWORD(&result, 0)&0x60000000)==0x60000000) { 3827 #elif QUAD 3828 if ((DFWORD(&result, 2)&0xffffff00)!=0 3829 || DFWORD(&result, 1)!=0 3830 || (DFWORD(&result, 0)&0x1c003fff)!=0 3831 || (DFWORD(&result, 0)&0x60000000)==0x60000000) { 3832 #endif 3833 set->status|=DEC_Invalid_operation; /* Invalid or out of range */ 3834 return 0; 3835 } 3836 /* get last twelve digits of the coefficent into hi & ho, base */ 3837 /* 10**9 (see GETCOEFFBILL): */ 3838 sourlo=DFWORD(&result, DECWORDS-1); 3839 lo=DPD2BIN0[sourlo&0x3ff] 3840 +DPD2BINK[(sourlo>>10)&0x3ff] 3841 +DPD2BINM[(sourlo>>20)&0x3ff]; 3842 sourpen=DFWORD(&result, DECWORDS-2); 3843 hi=DPD2BIN0[((sourpen<<2) | (sourlo>>30))&0x3ff]; 3844 3845 /* according to request, check range carefully */ 3846 if (unsign) { 3847 if (hi>4 || (hi==4 && lo>294967295) || (hi+lo!=0 && DFISSIGNED(&result))) { 3848 set->status|=DEC_Invalid_operation; /* out of range */ 3849 return 0; 3850 } 3851 return hi*BILLION+lo; 3852 } 3853 /* signed */ 3854 if (hi>2 || (hi==2 && lo>147483647)) { 3855 /* handle the usual edge case */ 3856 if (lo==147483648 && hi==2 && DFISSIGNED(&result)) return 0x80000000; 3857 set->status|=DEC_Invalid_operation; /* truly out of range */ 3858 return 0; 3859 } 3860 i=hi*BILLION+lo; 3861 if (DFISSIGNED(&result)) i=-i; 3862 return (uInt)i; 3863 } /* decToInt32 */ 3864 3865/* ------------------------------------------------------------------ */ 3866/* decToIntegral -- local routine to effect ToIntegral value */ 3867/* */ 3868/* result gets the result */ 3869/* df is the decFloat to round */ 3870/* set is the context */ 3871/* rmode is the rounding mode to use */ 3872/* exact is 1 if Inexact should be signalled */ 3873/* returns result */ 3874/* ------------------------------------------------------------------ */ 3875static decFloat * decToIntegral(decFloat *result, const decFloat *df, 3876 decContext *set, enum rounding rmode, 3877 Flag exact) { 3878 Int exp; /* exponent */ 3879 uInt sourhi; /* top word from source decFloat */ 3880 enum rounding saveround; /* saver */ 3881 uInt savestatus; /* .. */ 3882 decFloat zero; /* work */ 3883 3884 /* Start decoding the argument */ 3885 sourhi=DFWORD(df, 0); /* top word */ 3886 exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */ 3887 3888 if (EXPISSPECIAL(exp)) { /* is special? */ 3889 /* NaNs are handled as usual */ 3890 if (DFISNAN(df)) return decNaNs(result, df, NULL, set); 3891 /* must be infinite; return canonical infinity with sign of df */ 3892 return decInfinity(result, df); 3893 } 3894 3895 /* Here when the argument is finite */ 3896 /* complete extraction of the exponent */ 3897 exp+=GETECON(df)-DECBIAS; /* .. + continuation and unbias */ 3898 3899 if (exp>=0) return decCanonical(result, df); /* already integral */ 3900 3901 saveround=set->round; /* save rounding mode .. */ 3902 savestatus=set->status; /* .. and status */ 3903 set->round=rmode; /* set mode */ 3904 decFloatZero(&zero); /* make 0E+0 */ 3905 decFloatQuantize(result, df, &zero, set); /* 'integrate'; cannot fail */ 3906 set->round=saveround; /* restore rounding mode .. */ 3907 if (!exact) set->status=savestatus; /* .. and status, unless exact */ 3908 return result; 3909 } /* decToIntegral */ 3910