libdecnumber/decBasic.c

   1 /* Common base code for the decNumber C Library.
   2    Copyright (C) 2007, 2009 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
  17 Under Section 7 of GPL version 3, you are granted additional
  18 permissions described in the GCC Runtime Library Exception, version
  19 3.1, as published by the Free Software Foundation.
  20
  21 You should have received a copy of the GNU General Public License and
  22 a copy of the GCC Runtime Library Exception along with this program;
  23 see 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) */
  60 static decFloat *decDivide(decFloat *, const decFloat *,
  61                               const decFloat *, decContext *, uInt);
  62 static decFloat *decCanonical(decFloat *, const decFloat *);
  63 static void      decFiniteMultiply(bcdnum *, uByte *, const decFloat *,
  64                               const decFloat *);
  65 static decFloat *decInfinity(decFloat *, const decFloat *);
  66 static decFloat *decInvalid(decFloat *, decContext *);
  67 static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *,
  68                               decContext *);
  69 static Int       decNumCompare(const decFloat *, const decFloat *, Flag);
  70 static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *,
  71                               enum rounding, Flag);
  72 static 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 /* ------------------------------------------------------------------ */
  85 static 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) */
 177 static 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(&quotient, &num, set);         /* lay out the integer so far */
 666   DFWORD(&quotient, 0)^=DECFLOAT_Sign;       /* negate it */
 667   sign=DFWORD(dfl, 0);                       /* save sign of dfl */
 668   decFloatFMA(result, &quotient, dfr, dfl, set);
 669   if (!DFISZERO(result)) return result;
 670   /* if the result is zero the sign shall be sign of dfl */
 671   DFWORD(&quotient, 0)=sign;                 /* construct decFloat of sign */
 672   return decFloatCopySign(result, result, &quotient);
 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
 718 static 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 /* ------------------------------------------------------------------ */
1079 decFloat * 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). */
1102 const 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 */
1109 extern const Int DECTESTMSD[64];
1110 #endif
1111
1112 decFloat * 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 /* ------------------------------------------------------------------ */
1592 decFloat * 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 /* ------------------------------------------------------------------ */
1621 decFloat * 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 /* ------------------------------------------------------------------ */
1631 enum 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 /* ------------------------------------------------------------------ */
1663 const 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 /* ------------------------------------------------------------------ */
1687 decFloat * 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 /* ------------------------------------------------------------------ */
1711 decFloat * 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 /* ------------------------------------------------------------------ */
1737 decFloat * 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 /* ------------------------------------------------------------------ */
1803 decFloat * 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 /* ------------------------------------------------------------------ */
1827 decFloat * 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 /* ------------------------------------------------------------------ */
1841 decFloat * 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 /* ------------------------------------------------------------------ */
1856 decFloat * 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 /* ------------------------------------------------------------------ */
1872 decFloat * 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
1901 uInt 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. */
1973 decFloat * 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 /* ------------------------------------------------------------------ */
1989 decFloat * 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 /* ------------------------------------------------------------------ */
2006 decFloat * 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 /* ------------------------------------------------------------------ */
2284 decFloat * 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 /* ------------------------------------------------------------------ */
2319 decFloat * 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 /* ------------------------------------------------------------------ */
2349 decFloat * 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 /* ------------------------------------------------------------------ */
2377 uInt 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
2419 uInt decFloatIsFinite(const decFloat *df) {
2420   return !DFISSPECIAL(df);
2421   }
2422 uInt decFloatIsInfinite(const decFloat *df) {
2423   return DFISINF(df);
2424   }
2425 uInt decFloatIsInteger(const decFloat *df) {
2426   return DFISINT(df);
2427   }
2428 uInt decFloatIsNaN(const decFloat *df) {
2429   return DFISNAN(df);
2430   }
2431 uInt 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   }
2440 uInt decFloatIsSignaling(const decFloat *df) {
2441   return DFISSNAN(df);
2442   }
2443 uInt decFloatIsSignalling(const decFloat *df) {
2444   return DFISSNAN(df);
2445   }
2446 uInt decFloatIsSigned(const decFloat *df) {
2447   return DFISSIGNED(df);
2448   }
2449 uInt 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   }
2457 uInt 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 /* ------------------------------------------------------------------ */
2476 decFloat * 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 /* ------------------------------------------------------------------ */
2520 decFloat * 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 /* ------------------------------------------------------------------ */
2553 decFloat * 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 /* ------------------------------------------------------------------ */
2579 decFloat * 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 /* ------------------------------------------------------------------ */
2612 decFloat * 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 /* ------------------------------------------------------------------ */
2641 decFloat * 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 /* ------------------------------------------------------------------ */
2660 decFloat * 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 /* ------------------------------------------------------------------ */
2693 decFloat * 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 /* ------------------------------------------------------------------ */
2736 decFloat * 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 /* ------------------------------------------------------------------ */
2782 decFloat * 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 /* ------------------------------------------------------------------ */
2846 decFloat * 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 /* ------------------------------------------------------------------ */
2880 decFloat * 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 /* ------------------------------------------------------------------ */
2900 decFloat * 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;