libdecnumber/decBasic.c

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