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         ub = bufl;                           /* avoid type-pun violation */
1664         USHORTAT(ub)=0;
1665         uc = bufr;                           /* avoid type-pun violation */
1666         USHORTAT(uc)=0;
1667       #endif
1668       GETCOEFF(dfl, bufl+QUAD*2);            /* decode from decFloat */
1669       GETCOEFF(dfr, bufr+QUAD*2);            /* .. */
1670       /* all multiples of four, here */
1671       comp=0;                                /* assume equal */
1672       for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
1673         if (UINTAT(ub)==UINTAT(uc)) continue; /* so far so same */
1674         /* about to find a winner; go by bytes in case little-endian */
1675         for (;; ub++, uc++) {
1676           if (*ub==*uc) continue;
1677           if (*ub>*uc) comp=sigl;            /* difference found */
1678            else comp=-sigl;                  /* .. */
1679            break;
1680           }
1681         }
1682       } /* same NaN type and sign */
1683     }
1684    else {
1685     /* numeric comparison needed */
1686     comp=decNumCompare(dfl, dfr, 1);    /* total ordering */
1687     }
1688   decFloatZero(result);
1689   if (comp==0) return result;
1690   DFBYTE(result, DECBYTES-1)=0x01;      /* LSD=1 */
1691   if (comp<0) DFBYTE(result, 0)|=0x80;  /* set sign bit */
1692   return result;
1693   } /* decFloatCompareTotal */
1694
1695 /* ------------------------------------------------------------------ */
1696 /* decFloatCompareTotalMag -- compare magnitudes with total ordering  */
1697 /*                                                                    */
1698 /*   result gets the result of comparing abs(dfl) and abs(dfr)        */
1699 /*   dfl    is the first decFloat (lhs)                               */
1700 /*   dfr    is the second decFloat (rhs)                              */
1701 /*   returns result, which may be -1, 0, or 1                         */
1702 /* ------------------------------------------------------------------ */
1703 decFloat * decFloatCompareTotalMag(decFloat *result,
1704                                 const decFloat *dfl, const decFloat *dfr) {
1705   decFloat a, b;                        /* for copy if needed */
1706   /* copy and redirect signed operand(s) */
1707   if (DFISSIGNED(dfl)) {
1708     decFloatCopyAbs(&a, dfl);
1709     dfl=&a;
1710     }
1711   if (DFISSIGNED(dfr)) {
1712     decFloatCopyAbs(&b, dfr);
1713     dfr=&b;
1714     }
1715   return decFloatCompareTotal(result, dfl, dfr);
1716   } /* decFloatCompareTotalMag */
1717
1718 /* ------------------------------------------------------------------ */
1719 /* decFloatCopy -- copy a decFloat as-is                              */
1720 /*                                                                    */
1721 /*   result gets the copy of dfl                                      */
1722 /*   dfl    is the decFloat to copy                                   */
1723 /*   returns result                                                   */
1724 /*                                                                    */
1725 /* This is a bitwise operation; no errors or exceptions are possible. */
1726 /* ------------------------------------------------------------------ */
1727 decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) {
1728   if (dfl!=result) *result=*dfl;             /* copy needed */
1729   return result;
1730   } /* decFloatCopy */
1731
1732 /* ------------------------------------------------------------------ */
1733 /* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0     */
1734 /*                                                                    */
1735 /*   result gets the copy of dfl with sign bit 0                      */
1736 /*   dfl    is the decFloat to copy                                   */
1737 /*   returns result                                                   */
1738 /*                                                                    */
1739 /* This is a bitwise operation; no errors or exceptions are possible. */
1740 /* ------------------------------------------------------------------ */
1741 decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) {
1742   if (dfl!=result) *result=*dfl;        /* copy needed */
1743   DFBYTE(result, 0)&=~0x80;             /* zero sign bit */
1744   return result;
1745   } /* decFloatCopyAbs */
1746
1747 /* ------------------------------------------------------------------ */
1748 /* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
1749 /*                                                                    */
1750 /*   result gets the copy of dfl with sign bit inverted               */
1751 /*   dfl    is the decFloat to copy                                   */
1752 /*   returns result                                                   */
1753 /*                                                                    */
1754 /* This is a bitwise operation; no errors or exceptions are possible. */
1755 /* ------------------------------------------------------------------ */
1756 decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) {
1757   if (dfl!=result) *result=*dfl;        /* copy needed */
1758   DFBYTE(result, 0)^=0x80;              /* invert sign bit */
1759   return result;
1760   } /* decFloatCopyNegate */
1761
1762 /* ------------------------------------------------------------------ */
1763 /* decFloatCopySign -- copy a decFloat with the sign of another       */
1764 /*                                                                    */
1765 /*   result gets the result of copying dfl with the sign of dfr       */
1766 /*   dfl    is the first decFloat (lhs)                               */
1767 /*   dfr    is the second decFloat (rhs)                              */
1768 /*   returns result                                                   */
1769 /*                                                                    */
1770 /* This is a bitwise operation; no errors or exceptions are possible. */
1771 /* ------------------------------------------------------------------ */
1772 decFloat * decFloatCopySign(decFloat *result,
1773                             const decFloat *dfl, const decFloat *dfr) {
1774   uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80);   /* save sign bit */
1775   if (dfl!=result) *result=*dfl;             /* copy needed */
1776   DFBYTE(result, 0)&=~0x80;                  /* clear sign .. */
1777   DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* .. and set saved */
1778   return result;
1779   } /* decFloatCopySign */
1780
1781 /* ------------------------------------------------------------------ */
1782 /* decFloatDigits -- return the number of digits in a decFloat        */
1783 /*                                                                    */
1784 /*   df is the decFloat to investigate                                */
1785 /*   returns the number of significant digits in the decFloat; a      */
1786 /*     zero coefficient returns 1 as does an infinity (a NaN returns  */
1787 /*     the number of digits in the payload)                           */
1788 /* ------------------------------------------------------------------ */
1789 /* private macro to extract a declet according to provided formula */
1790 /* (form), and if it is non-zero then return the calculated digits */
1791 /* depending on the declet number (n), where n=0 for the most */
1792 /* significant declet; uses uInt dpd for work */
1793 #define dpdlenchk(n, form) {dpd=(form)&0x3ff;     \
1794   if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
1795 /* next one is used when it is known that the declet must be */
1796 /* non-zero, or is the final zero declet */
1797 #define dpdlendun(n, form) {dpd=(form)&0x3ff;     \
1798   if (dpd==0) return 1;                           \
1799   return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
1800
1801 uInt decFloatDigits(const decFloat *df) {
1802   uInt dpd;                        /* work */
1803   uInt sourhi=DFWORD(df, 0);       /* top word from source decFloat */
1804   #if QUAD
1805   uInt sourmh, sourml;
1806   #endif
1807   uInt sourlo;
1808
1809   if (DFISINF(df)) return 1;
1810   /* A NaN effectively has an MSD of 0; otherwise if non-zero MSD */
1811   /* then the coefficient is full-length */
1812   if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX;
1813
1814   #if DOUBLE
1815     if (sourhi&0x0003ffff) {     /* ends in first */
1816       dpdlenchk(0, sourhi>>8);
1817       sourlo=DFWORD(df, 1);
1818       dpdlendun(1, (sourhi<<2) | (sourlo>>30));
1819       } /* [cannot drop through] */
1820     sourlo=DFWORD(df, 1);  /* sourhi not involved now */
1821     if (sourlo&0xfff00000) {     /* in one of first two */
1822       dpdlenchk(1, sourlo>>30);  /* very rare */
1823       dpdlendun(2, sourlo>>20);
1824       } /* [cannot drop through] */
1825     dpdlenchk(3, sourlo>>10);
1826     dpdlendun(4, sourlo);
1827     /* [cannot drop through] */
1828
1829   #elif QUAD
1830     if (sourhi&0x00003fff) {     /* ends in first */
1831       dpdlenchk(0, sourhi>>4);
1832       sourmh=DFWORD(df, 1);
1833       dpdlendun(1, ((sourhi)<<6) | (sourmh>>26));
1834       } /* [cannot drop through] */
1835     sourmh=DFWORD(df, 1);
1836     if (sourmh) {
1837       dpdlenchk(1, sourmh>>26);
1838       dpdlenchk(2, sourmh>>16);
1839       dpdlenchk(3, sourmh>>6);
1840       sourml=DFWORD(df, 2);
1841       dpdlendun(4, ((sourmh)<<4) | (sourml>>28));
1842       } /* [cannot drop through] */
1843     sourml=DFWORD(df, 2);
1844     if (sourml) {
1845       dpdlenchk(4, sourml>>28);
1846       dpdlenchk(5, sourml>>18);
1847       dpdlenchk(6, sourml>>8);
1848       sourlo=DFWORD(df, 3);
1849       dpdlendun(7, ((sourml)<<2) | (sourlo>>30));
1850       } /* [cannot drop through] */
1851     sourlo=DFWORD(df, 3);
1852     if (sourlo&0xfff00000) {     /* in one of first two */
1853       dpdlenchk(7, sourlo>>30);  /* very rare */
1854       dpdlendun(8, sourlo>>20);
1855       } /* [cannot drop through] */
1856     dpdlenchk(9, sourlo>>10);
1857     dpdlendun(10, sourlo);
1858     /* [cannot drop through] */
1859   #endif
1860   } /* decFloatDigits */
1861
1862 /* ------------------------------------------------------------------ */
1863 /* decFloatDivide -- divide a decFloat by another                     */
1864 /*                                                                    */
1865 /*   result gets the result of dividing dfl by dfr:                   */
1866 /*   dfl    is the first decFloat (lhs)                               */
1867 /*   dfr    is the second decFloat (rhs)                              */
1868 /*   set    is the context                                            */
1869 /*   returns result                                                   */
1870 /*                                                                    */
1871 /* ------------------------------------------------------------------ */
1872 /* This is just a wrapper. */
1873 decFloat * decFloatDivide(decFloat *result,
1874                           const decFloat *dfl, const decFloat *dfr,
1875                           decContext *set) {
1876   return decDivide(result, dfl, dfr, set, DIVIDE);
1877   } /* decFloatDivide */
1878
1879 /* ------------------------------------------------------------------ */
1880 /* decFloatDivideInteger -- integer divide a decFloat by another      */
1881 /*                                                                    */
1882 /*   result gets the result of dividing dfl by dfr:                   */
1883 /*   dfl    is the first decFloat (lhs)                               */
1884 /*   dfr    is the second decFloat (rhs)                              */
1885 /*   set    is the context                                            */
1886 /*   returns result                                                   */
1887 /*                                                                    */
1888 /* ------------------------------------------------------------------ */
1889 decFloat * decFloatDivideInteger(decFloat *result,
1890                              const decFloat *dfl, const decFloat *dfr,
1891                              decContext *set) {
1892   return decDivide(result, dfl, dfr, set, DIVIDEINT);
1893   } /* decFloatDivideInteger */
1894
1895 /* ------------------------------------------------------------------ */
1896 /* decFloatFMA -- multiply and add three decFloats, fused             */
1897 /*                                                                    */
1898 /*   result gets the result of (dfl*dfr)+dff with a single rounding   */
1899 /*   dfl    is the first decFloat (lhs)                               */
1900 /*   dfr    is the second decFloat (rhs)                              */
1901 /*   dff    is the final decFloat (fhs)                               */
1902 /*   set    is the context                                            */
1903 /*   returns result                                                   */
1904 /*                                                                    */
1905 /* ------------------------------------------------------------------ */
1906 decFloat * decFloatFMA(decFloat *result, const decFloat *dfl,
1907                        const decFloat *dfr, const decFloat *dff,
1908                        decContext *set) {
1909   /* The accumulator has the bytes needed for FiniteMultiply, plus */
1910   /* one byte to the left in case of carry, plus DECPMAX+2 to the */
1911   /* right for the final addition (up to full fhs + round & sticky) */
1912   #define FMALEN (1+ (DECPMAX9*18) +DECPMAX+2)
1913   uByte  acc[FMALEN];              /* for multiplied coefficient in BCD */
1914                                    /* .. and for final result */
1915   bcdnum mul;                      /* for multiplication result */
1916   bcdnum fin;                      /* for final operand, expanded */
1917   uByte  coe[DECPMAX];             /* dff coefficient in BCD */
1918   bcdnum *hi, *lo;                 /* bcdnum with higher/lower exponent */
1919   uInt   diffsign;                 /* non-zero if signs differ */
1920   uInt   hipad;                    /* pad digit for hi if needed */
1921   Int    padding;                  /* excess exponent */
1922   uInt   carry;                    /* +1 for ten's complement and during add */
1923   uByte  *ub, *uh, *ul;            /* work */
1924
1925   /* handle all the special values [any special operand leads to a */
1926   /* special result] */
1927   if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr) || DFISSPECIAL(dff)) {
1928     decFloat proxy;                /* multiplication result proxy */
1929     /* NaNs are handled as usual, giving priority to sNaNs */
1930     if (DFISSNAN(dfl) || DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
1931     if (DFISSNAN(dff)) return decNaNs(result, dff, NULL, set);
1932     if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
1933     if (DFISNAN(dff)) return decNaNs(result, dff, NULL, set);
1934     /* One or more of the three is infinite */
1935     /* infinity times zero is bad */
1936     decFloatZero(&proxy);
1937     if (DFISINF(dfl)) {
1938       if (DFISZERO(dfr)) return decInvalid(result, set);
1939       decInfinity(&proxy, &proxy);
1940       }
1941      else if (DFISINF(dfr)) {
1942       if (DFISZERO(dfl)) return decInvalid(result, set);
1943       decInfinity(&proxy, &proxy);
1944       }
1945     /* compute sign of multiplication and place in proxy */
1946     DFWORD(&proxy, 0)|=(DFWORD(dfl, 0)^DFWORD(dfr, 0))&DECFLOAT_Sign;
1947     if (!DFISINF(dff)) return decFloatCopy(result, &proxy);
1948     /* dff is Infinite */
1949     if (!DFISINF(&proxy)) return decInfinity(result, dff);
1950     /* both sides of addition are infinite; different sign is bad */
1951     if ((DFWORD(dff, 0)&DECFLOAT_Sign)!=(DFWORD(&proxy, 0)&DECFLOAT_Sign))
1952       return decInvalid(result, set);
1953     return decFloatCopy(result, &proxy);
1954     }
1955
1956   /* Here when all operands are finite */
1957
1958   /* First multiply dfl*dfr */
1959   decFiniteMultiply(&mul, acc+1, dfl, dfr);
1960   /* The multiply is complete, exact and unbounded, and described in */
1961   /* mul with the coefficient held in acc[1...] */
1962
1963   /* now add in dff; the algorithm is essentially the same as */
1964   /* decFloatAdd, but the code is different because the code there */
1965   /* is highly optimized for adding two numbers of the same size */
1966   fin.exponent=GETEXPUN(dff);           /* get dff exponent and sign */
1967   fin.sign=DFWORD(dff, 0)&DECFLOAT_Sign;
1968   diffsign=mul.sign^fin.sign;           /* note if signs differ */
1969   fin.msd=coe;
1970   fin.lsd=coe+DECPMAX-1;
1971   GETCOEFF(dff, coe);                   /* extract the coefficient */
1972
1973   /* now set hi and lo so that hi points to whichever of mul and fin */
1974   /* has the higher exponent and lo point to the other [don't care if */
1975   /* the same] */
1976   if (mul.exponent>=fin.exponent) {
1977     hi=&mul;
1978     lo=&fin;
1979     }
1980    else {
1981     hi=&fin;
1982     lo=&mul;
1983     }
1984
1985   /* remove leading zeros on both operands; this will save time later */
1986   /* and make testing for zero trivial */
1987   for (; UINTAT(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4;
1988   for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++;
1989   for (; UINTAT(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
1990   for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
1991
1992   /* if hi is zero then result will be lo (which has the smaller */
1993   /* exponent), which also may need to be tested for zero for the */
1994   /* weird IEEE 754 sign rules */
1995   if (*hi->msd==0 && hi->msd==hi->lsd) {     /* hi is zero */
1996     /* "When the sum of two operands with opposite signs is */
1997     /* exactly zero, the sign of that sum shall be '+' in all */
1998     /* rounding modes except round toward -Infinity, in which */
1999     /* mode that sign shall be '-'." */
2000     if (diffsign) {
2001       if (*lo->msd==0 && lo->msd==lo->lsd) { /* lo is zero */
2002         lo->sign=0;
2003         if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
2004         } /* diffsign && lo=0 */
2005       } /* diffsign */
2006     return decFinalize(result, lo, set);     /* may need clamping */
2007     } /* numfl is zero */
2008   /* [here, both are minimal length and hi is non-zero] */
2009
2010   /* if signs differ, take the ten's complement of hi (zeros to the */
2011   /* right do not matter because the complement of zero is zero); */
2012   /* the +1 is done later, as part of the addition, inserted at the */
2013   /* correct digit */
2014   hipad=0;
2015   carry=0;
2016   if (diffsign) {
2017     hipad=9;
2018     carry=1;
2019     /* exactly the correct number of digits must be inverted */
2020     for (uh=hi->msd; uh<hi->lsd-3; uh+=4) UINTAT(uh)=0x09090909-UINTAT(uh);
2021     for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh);
2022     }
2023
2024   /* ready to add; note that hi has no leading zeros so gap */
2025   /* calculation does not have to be as pessimistic as in decFloatAdd */
2026   /* (this is much more like the arbitrary-precision algorithm in */
2027   /* Rexx and decNumber) */
2028
2029   /* padding is the number of zeros that would need to be added to hi */
2030   /* for its lsd to be aligned with the lsd of lo */
2031   padding=hi->exponent-lo->exponent;
2032   /* printf("FMA pad %ld\n", (LI)padding); */
2033
2034   /* the result of the addition will be built into the accumulator, */
2035   /* starting from the far right; this could be either hi or lo */
2036   ub=acc+FMALEN-1;                 /* where lsd of result will go */
2037   ul=lo->lsd;                      /* lsd of rhs */
2038
2039   if (padding!=0) {                /* unaligned */
2040     /* if the msd of lo is more than DECPMAX+2 digits to the right of */
2041     /* the original msd of hi then it can be reduced to a single */
2042     /* digit at the right place, as it stays clear of hi digits */
2043     /* [it must be DECPMAX+2 because during a subtraction the msd */
2044     /* could become 0 after a borrow from 1.000 to 0.9999...] */
2045     Int hilen=(Int)(hi->lsd-hi->msd+1); /* lengths */
2046     Int lolen=(Int)(lo->lsd-lo->msd+1); /* .. */
2047     Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
2048     Int reduce=newexp-lo->exponent;
2049     if (reduce>0) {                     /* [= case gives reduce=0 nop] */
2050       /* printf("FMA reduce: %ld\n", (LI)reduce); */
2051       if (reduce>=lolen) {              /* eating all */
2052         lo->lsd=lo->msd;                /* reduce to single digit */
2053         lo->exponent=newexp;            /* [known to be non-zero] */
2054         }
2055        else { /* < */
2056         uByte *up=lo->lsd;
2057         lo->lsd=lo->lsd-reduce;
2058         if (*lo->lsd==0)                /* could need sticky bit */
2059          for (; up>lo->lsd; up--) {     /* search discarded digits */
2060           if (*up!=0) {                 /* found one... */
2061             *lo->lsd=1;                 /* set sticky bit */
2062             break;
2063             }
2064           }
2065         lo->exponent+=reduce;
2066         }
2067       padding=hi->exponent-lo->exponent; /* recalculate */
2068       ul=lo->lsd;                        /* .. */
2069       } /* maybe reduce */
2070     /* padding is now <= DECPMAX+2 but still > 0; tricky DOUBLE case */
2071     /* is when hi is a 1 that will become a 0.9999... by subtraction: */
2072     /*   hi:   1                                   E+16 */
2073     /*   lo:    .................1000000000000000  E-16 */
2074     /* which for the addition pads and reduces to: */
2075     /*   hi:   1000000000000000000                 E-2 */
2076     /*   lo:    .................1                 E-2 */
2077     #if DECCHECK
2078       if (padding>DECPMAX+2) printf("FMA excess padding: %ld\n", (LI)padding);
2079       if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding);
2080       /* printf("FMA padding: %ld\n", (LI)padding); */
2081     #endif
2082     /* padding digits can now be set in the result; one or more of */
2083     /* these will come from lo; others will be zeros in the gap */
2084     for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul;
2085     for (;padding>0; padding--, ub--) *ub=0; /* mind the gap */
2086     }
2087
2088   /* addition now complete to the right of the rightmost digit of hi */
2089   uh=hi->lsd;
2090
2091   /* carry was set up depending on ten's complement above; do the add... */
2092   for (;; ub--) {
2093     uInt hid, lod;
2094     if (uh<hi->msd) {
2095       if (ul<lo->msd) break;
2096       hid=hipad;
2097       }
2098      else hid=*uh--;
2099     if (ul<lo->msd) lod=0;
2100      else lod=*ul--;
2101     *ub=(uByte)(carry+hid+lod);
2102     if (*ub<10) carry=0;
2103      else {
2104       *ub-=10;
2105       carry=1;
2106       }
2107     } /* addition loop */
2108
2109   /* addition complete -- now handle carry, borrow, etc. */
2110   /* use lo to set up the num (its exponent is already correct, and */
2111   /* sign usually is) */
2112   lo->msd=ub+1;
2113   lo->lsd=acc+FMALEN-1;
2114   /* decShowNum(lo, "lo"); */
2115   if (!diffsign) {                 /* same-sign addition */
2116     if (carry) {                   /* carry out */
2117       *ub=1;                       /* place the 1 .. */
2118       lo->msd--;                   /* .. and update */
2119       }
2120     } /* same sign */
2121    else {                          /* signs differed (subtraction) */
2122     if (!carry) {                  /* no carry out means hi<lo */
2123       /* borrowed -- take ten's complement of the right digits */
2124       lo->sign=hi->sign;           /* sign is lhs sign */
2125       for (ul=lo->msd; ul<lo->lsd-3; ul+=4) UINTAT(ul)=0x09090909-UINTAT(ul);
2126       for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); /* [leaves ul at lsd+1] */
2127       /* complete the ten's complement by adding 1 [cannot overrun] */
2128       for (ul--; *ul==9; ul--) *ul=0;
2129       *ul+=1;
2130       } /* borrowed */
2131      else {                        /* carry out means hi>=lo */
2132       /* sign to use is lo->sign */
2133       /* all done except for the special IEEE 754 exact-zero-result */
2134       /* rule (see above); while testing for zero, strip leading */
2135       /* zeros (which will save decFinalize doing it) */
2136       for (; UINTAT(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
2137       for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
2138       if (*lo->msd==0) {           /* must be true zero (and diffsign) */
2139         lo->sign=0;                /* assume + */
2140         if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
2141         }
2142       /* [else was not zero, might still have leading zeros] */
2143       } /* subtraction gave positive result */
2144     } /* diffsign */
2145
2146   return decFinalize(result, lo, set);  /* round, check, and lay out */
2147   } /* decFloatFMA */
2148
2149 /* ------------------------------------------------------------------ */
2150 /* decFloatFromInt -- initialise a decFloat from an Int               */
2151 /*                                                                    */
2152 /*   result gets the converted Int                                    */
2153 /*   n      is the Int to convert                                     */
2154 /*   returns result                                                   */
2155 /*                                                                    */
2156 /* The result is Exact; no errors or exceptions are possible.         */
2157 /* ------------------------------------------------------------------ */
2158 decFloat * decFloatFromInt32(decFloat *result, Int n) {
2159   uInt u=(uInt)n;                       /* copy as bits */
2160   uInt encode;                          /* work */
2161   DFWORD(result, 0)=ZEROWORD;           /* always */
2162   #if QUAD
2163     DFWORD(result, 1)=0;
2164     DFWORD(result, 2)=0;
2165   #endif
2166   if (n<0) {                            /* handle -n with care */
2167     /* [This can be done without the test, but is then slightly slower] */
2168     u=(~u)+1;
2169     DFWORD(result, 0)|=DECFLOAT_Sign;
2170     }
2171   /* Since the maximum value of u now is 2**31, only the low word of */
2172   /* result is affected */
2173   encode=BIN2DPD[u%1000];
2174   u/=1000;
2175   encode|=BIN2DPD[u%1000]<<10;
2176   u/=1000;
2177   encode|=BIN2DPD[u%1000]<<20;
2178   u/=1000;                              /* now 0, 1, or 2 */
2179   encode|=u<<30;
2180   DFWORD(result, DECWORDS-1)=encode;
2181   return result;
2182   } /* decFloatFromInt32 */
2183
2184 /* ------------------------------------------------------------------ */
2185 /* decFloatFromUInt -- initialise a decFloat from a uInt              */
2186 /*                                                                    */
2187 /*   result gets the converted uInt                                   */
2188 /*   n      is the uInt to convert                                    */
2189 /*   returns result                                                   */
2190 /*                                                                    */
2191 /* The result is Exact; no errors or exceptions are possible.         */
2192 /* ------------------------------------------------------------------ */
2193 decFloat * decFloatFromUInt32(decFloat *result, uInt u) {
2194   uInt encode;                          /* work */
2195   DFWORD(result, 0)=ZEROWORD;           /* always */
2196   #if QUAD
2197     DFWORD(result, 1)=0;
2198     DFWORD(result, 2)=0;
2199   #endif
2200   encode=BIN2DPD[u%1000];
2201   u/=1000;
2202   encode|=BIN2DPD[u%1000]<<10;
2203   u/=1000;
2204   encode|=BIN2DPD[u%1000]<<20;
2205   u/=1000;                              /* now 0 -> 4 */
2206   encode|=u<<30;
2207   DFWORD(result, DECWORDS-1)=encode;
2208   DFWORD(result, DECWORDS-2)|=u>>2;     /* rarely non-zero */
2209   return result;
2210   } /* decFloatFromUInt32 */
2211
2212 /* ------------------------------------------------------------------ */
2213 /* decFloatInvert -- logical digitwise INVERT of a decFloat           */
2214 /*                                                                    */
2215 /*   result gets the result of INVERTing df                           */
2216 /*   df     is the decFloat to invert                                 */
2217 /*   set    is the context                                            */
2218 /*   returns result, which will be canonical with sign=0              */
2219 /*                                                                    */
2220 /* The operand must be positive, finite with exponent q=0, and        */
2221 /* comprise just zeros and ones; if not, Invalid operation results.   */
2222 /* ------------------------------------------------------------------ */
2223 decFloat * decFloatInvert(decFloat *result, const decFloat *df,
2224                           decContext *set) {
2225   uInt sourhi=DFWORD(df, 0);            /* top word of dfs */
2226
2227   if (!DFISUINT01(df) || !DFISCC01(df)) return decInvalid(result, set);
2228   /* the operand is a finite integer (q=0) */
2229   #if DOUBLE
2230    DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04009124);
2231    DFWORD(result, 1)=(~DFWORD(df, 1))   &0x49124491;
2232   #elif QUAD
2233    DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04000912);
2234    DFWORD(result, 1)=(~DFWORD(df, 1))   &0x44912449;
2235    DFWORD(result, 2)=(~DFWORD(df, 2))   &0x12449124;
2236    DFWORD(result, 3)=(~DFWORD(df, 3))   &0x49124491;
2237   #endif
2238   return result;
2239   } /* decFloatInvert */
2240
2241 /* ------------------------------------------------------------------ */
2242 /* decFloatIs -- decFloat tests (IsSigned, etc.)                      */
2243 /*                                                                    */
2244 /*   df is the decFloat to test                                       */
2245 /*   returns 0 or 1 in an int32_t                                     */
2246 /*                                                                    */
2247 /* Many of these could be macros, but having them as real functions   */
2248 /* is a bit cleaner (and they can be referred to here by the generic  */
2249 /* names)                                                             */
2250 /* ------------------------------------------------------------------ */
2251 uInt decFloatIsCanonical(const decFloat *df) {
2252   if (DFISSPECIAL(df)) {
2253     if (DFISINF(df)) {
2254       if (DFWORD(df, 0)&ECONMASK) return 0;  /* exponent continuation */
2255       if (!DFISCCZERO(df)) return 0;         /* coefficient continuation */
2256       return 1;
2257       }
2258     /* is a NaN */
2259     if (DFWORD(df, 0)&ECONNANMASK) return 0; /* exponent continuation */
2260     if (DFISCCZERO(df)) return 1;            /* coefficient continuation */
2261     /* drop through to check payload */
2262     }
2263   { /* declare block */
2264   #if DOUBLE
2265     uInt sourhi=DFWORD(df, 0);
2266     uInt sourlo=DFWORD(df, 1);
2267     if (CANONDPDOFF(sourhi, 8)
2268      && CANONDPDTWO(sourhi, sourlo, 30)
2269      && CANONDPDOFF(sourlo, 20)
2270      && CANONDPDOFF(sourlo, 10)
2271      && CANONDPDOFF(sourlo, 0)) return 1;
2272   #elif QUAD
2273     uInt sourhi=DFWORD(df, 0);
2274     uInt sourmh=DFWORD(df, 1);
2275     uInt sourml=DFWORD(df, 2);
2276     uInt sourlo=DFWORD(df, 3);
2277     if (CANONDPDOFF(sourhi, 4)
2278      && CANONDPDTWO(sourhi, sourmh, 26)
2279      && CANONDPDOFF(sourmh, 16)
2280      && CANONDPDOFF(sourmh, 6)
2281      && CANONDPDTWO(sourmh, sourml, 28)
2282      && CANONDPDOFF(sourml, 18)
2283      && CANONDPDOFF(sourml, 8)
2284      && CANONDPDTWO(sourml, sourlo, 30)
2285      && CANONDPDOFF(sourlo, 20)
2286      && CANONDPDOFF(sourlo, 10)
2287      && CANONDPDOFF(sourlo, 0)) return 1;
2288   #endif
2289   } /* block */
2290   return 0;    /* a declet is non-canonical */
2291   }
2292
2293 uInt decFloatIsFinite(const decFloat *df) {
2294   return !DFISSPECIAL(df);
2295   }
2296 uInt decFloatIsInfinite(const decFloat *df) {
2297   return DFISINF(df);
2298   }
2299 uInt decFloatIsInteger(const decFloat *df) {
2300   return DFISINT(df);
2301   }
2302 uInt decFloatIsNaN(const decFloat *df) {
2303   return DFISNAN(df);
2304   }
2305 uInt decFloatIsNormal(const decFloat *df) {
2306   Int exp;                         /* exponent */
2307   if (DFISSPECIAL(df)) return 0;
2308   if (DFISZERO(df)) return 0;
2309   /* is finite and non-zero */
2310   exp=GETEXPUN(df)                 /* get unbiased exponent .. */
2311      +decFloatDigits(df)-1;        /* .. and make adjusted exponent */
2312   return (exp>=DECEMIN);           /* < DECEMIN is subnormal */
2313   }
2314 uInt decFloatIsSignaling(const decFloat *df) {
2315   return DFISSNAN(df);
2316   }
2317 uInt decFloatIsSignalling(const decFloat *df) {
2318   return DFISSNAN(df);
2319   }
2320 uInt decFloatIsSigned(const decFloat *df) {
2321   return DFISSIGNED(df);
2322   }
2323 uInt decFloatIsSubnormal(const decFloat *df) {
2324   if (DFISSPECIAL(df)) return 0;
2325   /* is finite */
2326   if (decFloatIsNormal(df)) return 0;
2327   /* it is <Nmin, but could be zero */
2328   if (DFISZERO(df)) return 0;
2329   return 1;                                  /* is subnormal */
2330   }
2331 uInt decFloatIsZero(const decFloat *df) {
2332   return DFISZERO(df);
2333   } /* decFloatIs... */
2334
2335 /* ------------------------------------------------------------------ */
2336 /* decFloatLogB -- return adjusted exponent, by 754r rules            */
2337 /*                                                                    */
2338 /*   result gets the adjusted exponent as an integer, or a NaN etc.   */
2339 /*   df     is the decFloat to be examined                            */
2340 /*   set    is the context                                            */
2341 /*   returns result                                                   */
2342 /*                                                                    */
2343 /* Notable cases:                                                     */
2344 /*   A<0 -> Use |A|                                                   */
2345 /*   A=0 -> -Infinity (Division by zero)                              */
2346 /*   A=Infinite -> +Infinity (Exact)                                  */
2347 /*   A=1 exactly -> 0 (Exact)                                         */
2348 /*   NaNs are propagated as usual                                     */
2349 /* ------------------------------------------------------------------ */
2350 decFloat * decFloatLogB(decFloat *result, const decFloat *df,
2351                         decContext *set) {
2352   Int ae;                                    /* adjusted exponent */
2353   if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2354   if (DFISINF(df)) {
2355     DFWORD(result, 0)=0;                     /* need +ve */
2356     return decInfinity(result, result);      /* canonical +Infinity */
2357     }
2358   if (DFISZERO(df)) {
2359     set->status|=DEC_Division_by_zero;       /* as per 754r */
2360     DFWORD(result, 0)=DECFLOAT_Sign;         /* make negative */
2361     return decInfinity(result, result);      /* canonical -Infinity */
2362     }
2363   ae=GETEXPUN(df)                       /* get unbiased exponent .. */
2364     +decFloatDigits(df)-1;              /* .. and make adjusted exponent */
2365   /* ae has limited range (3 digits for DOUBLE and 4 for QUAD), so */
2366   /* it is worth using a special case of decFloatFromInt32 */
2367   DFWORD(result, 0)=ZEROWORD;           /* always */
2368   if (ae<0) {
2369     DFWORD(result, 0)|=DECFLOAT_Sign;   /* -0 so far */
2370     ae=-ae;
2371     }
2372   #if DOUBLE
2373     DFWORD(result, 1)=BIN2DPD[ae];      /* a single declet */
2374   #elif QUAD
2375     DFWORD(result, 1)=0;
2376     DFWORD(result, 2)=0;
2377     DFWORD(result, 3)=(ae/1000)<<10;    /* is <10, so need no DPD encode */
2378     DFWORD(result, 3)|=BIN2DPD[ae%1000];
2379   #endif
2380   return result;
2381   } /* decFloatLogB */
2382
2383 /* ------------------------------------------------------------------ */
2384 /* decFloatMax -- return maxnum of two operands                       */
2385 /*                                                                    */
2386 /*   result gets the chosen decFloat                                  */
2387 /*   dfl    is the first decFloat (lhs)                               */
2388 /*   dfr    is the second decFloat (rhs)                              */
2389 /*   set    is the context                                            */
2390 /*   returns result                                                   */
2391 /*                                                                    */
2392 /* If just one operand is a quiet NaN it is ignored.                  */
2393 /* ------------------------------------------------------------------ */
2394 decFloat * decFloatMax(decFloat *result,
2395                        const decFloat *dfl, const decFloat *dfr,
2396                        decContext *set) {
2397   Int comp;
2398   if (DFISNAN(dfl)) {
2399     /* sNaN or both NaNs leads to normal NaN processing */
2400     if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
2401     return decCanonical(result, dfr);        /* RHS is numeric */
2402     }
2403   if (DFISNAN(dfr)) {
2404     /* sNaN leads to normal NaN processing (both NaN handled above) */
2405     if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2406     return decCanonical(result, dfl);        /* LHS is numeric */
2407     }
2408   /* Both operands are numeric; numeric comparison needed -- use */
2409   /* total order for a well-defined choice (and +0 > -0) */
2410   comp=decNumCompare(dfl, dfr, 1);
2411   if (comp>=0) return decCanonical(result, dfl);
2412   return decCanonical(result, dfr);
2413   } /* decFloatMax */
2414
2415 /* ------------------------------------------------------------------ */
2416 /* decFloatMaxMag -- return maxnummag of two operands                 */
2417 /*                                                                    */
2418 /*   result gets the chosen decFloat                                  */
2419 /*   dfl    is the first decFloat (lhs)                               */
2420 /*   dfr    is the second decFloat (rhs)                              */
2421 /*   set    is the context                                            */
2422 /*   returns result                                                   */
2423 /*                                                                    */
2424 /* Returns according to the magnitude comparisons if both numeric and */
2425 /* unequal, otherwise returns maxnum                                  */
2426 /* ------------------------------------------------------------------ */
2427 decFloat * decFloatMaxMag(decFloat *result,
2428                        const decFloat *dfl, const decFloat *dfr,
2429                        decContext *set) {
2430   Int comp;
2431   decFloat absl, absr;
2432   if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMax(result, dfl, dfr, set);
2433
2434   decFloatCopyAbs(&absl, dfl);
2435   decFloatCopyAbs(&absr, dfr);
2436   comp=decNumCompare(&absl, &absr, 0);
2437   if (comp>0) return decCanonical(result, dfl);
2438   if (comp<0) return decCanonical(result, dfr);
2439   return decFloatMax(result, dfl, dfr, set);
2440   } /* decFloatMaxMag */
2441
2442 /* ------------------------------------------------------------------ */
2443 /* decFloatMin -- return minnum of two operands                       */
2444 /*                                                                    */
2445 /*   result gets the chosen decFloat                                  */
2446 /*   dfl    is the first decFloat (lhs)                               */
2447 /*   dfr    is the second decFloat (rhs)                              */
2448 /*   set    is the context                                            */
2449 /*   returns result                                                   */
2450 /*                                                                    */
2451 /* If just one operand is a quiet NaN it is ignored.                  */
2452 /* ------------------------------------------------------------------ */
2453 decFloat * decFloatMin(decFloat *result,
2454                        const decFloat *dfl, const decFloat *dfr,
2455                        decContext *set) {
2456   Int comp;
2457   if (DFISNAN(dfl)) {
2458     /* sNaN or both NaNs leads to normal NaN processing */
2459     if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
2460     return decCanonical(result, dfr);        /* RHS is numeric */
2461     }
2462   if (DFISNAN(dfr)) {
2463     /* sNaN leads to normal NaN processing (both NaN handled above) */
2464     if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2465     return decCanonical(result, dfl);        /* LHS is numeric */
2466     }
2467   /* Both operands are numeric; numeric comparison needed -- use */
2468   /* total order for a well-defined choice (and +0 > -0) */
2469   comp=decNumCompare(dfl, dfr, 1);
2470   if (comp<=0) return decCanonical(result, dfl);
2471   return decCanonical(result, dfr);
2472   } /* decFloatMin */
2473
2474 /* ------------------------------------------------------------------ */
2475 /* decFloatMinMag -- return minnummag of two operands                 */
2476 /*                                                                    */
2477 /*   result gets the chosen decFloat                                  */
2478 /*   dfl    is the first decFloat (lhs)                               */
2479 /*   dfr    is the second decFloat (rhs)                              */
2480 /*   set    is the context                                            */
2481 /*   returns result                                                   */
2482 /*                                                                    */
2483 /* Returns according to the magnitude comparisons if both numeric and */
2484 /* unequal, otherwise returns minnum                                  */
2485 /* ------------------------------------------------------------------ */
2486 decFloat * decFloatMinMag(decFloat *result,
2487                        const decFloat *dfl, const decFloat *dfr,
2488                        decContext *set) {
2489   Int comp;
2490   decFloat absl, absr;
2491   if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMin(result, dfl, dfr, set);
2492
2493   decFloatCopyAbs(&absl, dfl);
2494   decFloatCopyAbs(&absr, dfr);
2495   comp=decNumCompare(&absl, &absr, 0);
2496   if (comp<0) return decCanonical(result, dfl);
2497   if (comp>0) return decCanonical(result, dfr);
2498   return decFloatMin(result, dfl, dfr, set);
2499   } /* decFloatMinMag */
2500
2501 /* ------------------------------------------------------------------ */
2502 /* decFloatMinus -- negate value, heeding NaNs, etc.                  */
2503 /*                                                                    */
2504 /*   result gets the canonicalized 0-df                               */
2505 /*   df     is the decFloat to minus                                  */
2506 /*   set    is the context                                            */
2507 /*   returns result                                                   */
2508 /*                                                                    */
2509 /* This has the same effect as 0-df where the exponent of the zero is */
2510 /* the same as that of df (if df is finite).                          */
2511 /* The effect is also the same as decFloatCopyNegate except that NaNs */
2512 /* are handled normally (the sign of a NaN is not affected, and an    */
2513 /* sNaN will signal), the result is canonical, and zero gets sign 0.  */
2514 /* ------------------------------------------------------------------ */
2515 decFloat * decFloatMinus(decFloat *result, const decFloat *df,
2516                          decContext *set) {
2517   if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2518   decCanonical(result, df);                       /* copy and check */
2519   if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80;     /* turn off sign bit */
2520    else DFBYTE(result, 0)^=0x80;                  /* flip sign bit */
2521   return result;
2522   } /* decFloatMinus */
2523
2524 /* ------------------------------------------------------------------ */
2525 /* decFloatMultiply -- multiply two decFloats                         */
2526 /*                                                                    */
2527 /*   result gets the result of multiplying dfl and dfr:               */
2528 /*   dfl    is the first decFloat (lhs)                               */
2529 /*   dfr    is the second decFloat (rhs)                              */
2530 /*   set    is the context                                            */
2531 /*   returns result                                                   */
2532 /*                                                                    */
2533 /* ------------------------------------------------------------------ */
2534 decFloat * decFloatMultiply(decFloat *result,
2535                             const decFloat *dfl, const decFloat *dfr,
2536                             decContext *set) {
2537   bcdnum num;                      /* for final conversion */
2538   uByte  bcdacc[DECPMAX9*18+1];    /* for coefficent in BCD */
2539
2540   if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
2541     /* NaNs are handled as usual */
2542     if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2543     /* infinity times zero is bad */
2544     if (DFISINF(dfl) && DFISZERO(dfr)) return decInvalid(result, set);
2545     if (DFISINF(dfr) && DFISZERO(dfl)) return decInvalid(result, set);
2546     /* both infinite; return canonical infinity with computed sign */
2547     DFWORD(result, 0)=DFWORD(dfl, 0)^DFWORD(dfr, 0); /* compute sign */
2548     return decInfinity(result, result);
2549     }
2550
2551   /* Here when both operands are finite */
2552   decFiniteMultiply(&num, bcdacc, dfl, dfr);
2553   return decFinalize(result, &num, set); /* round, check, and lay out */
2554   } /* decFloatMultiply */
2555
2556 /* ------------------------------------------------------------------ */
2557 /* decFloatNextMinus -- next towards -Infinity                        */
2558 /*                                                                    */
2559 /*   result gets the next lesser decFloat                             */
2560 /*   dfl    is the decFloat to start with                             */
2561 /*   set    is the context                                            */
2562 /*   returns result                                                   */
2563 /*                                                                    */
2564 /* This is 754r nextdown; Invalid is the only status possible (from   */
2565 /* an sNaN).                                                          */
2566 /* ------------------------------------------------------------------ */
2567 decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
2568                              decContext *set) {
2569   decFloat delta;                       /* tiny increment */
2570   uInt savestat;                        /* saves status */
2571   enum rounding saveround;              /* .. and mode */
2572
2573   /* +Infinity is the special case */
2574   if (DFISINF(dfl) && !DFISSIGNED(dfl)) {
2575     DFSETNMAX(result);
2576     return result;                      /* [no status to set] */
2577     }
2578   /* other cases are effected by sutracting a tiny delta -- this */
2579   /* should be done in a wider format as the delta is unrepresentable */
2580   /* here (but can be done with normal add if the sign of zero is */
2581   /* treated carefully, because no Inexactitude is interesting); */
2582   /* rounding to -Infinity then pushes the result to next below */
2583   decFloatZero(&delta);                 /* set up tiny delta */
2584   DFWORD(&delta, DECWORDS-1)=1;         /* coefficient=1 */
2585   DFWORD(&delta, 0)=DECFLOAT_Sign;      /* Sign=1 + biased exponent=0 */
2586   /* set up for the directional round */
2587   saveround=set->round;                 /* save mode */
2588   set->round=DEC_ROUND_FLOOR;           /* .. round towards -Infinity */
2589   savestat=set->status;                 /* save status */
2590   decFloatAdd(result, dfl, &delta, set);
2591   /* Add rules mess up the sign when going from +Ntiny to 0 */
2592   if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
2593   set->status&=DEC_Invalid_operation;   /* preserve only sNaN status */
2594   set->status|=savestat;                /* restore pending flags */
2595   set->round=saveround;                 /* .. and mode */
2596   return result;
2597   } /* decFloatNextMinus */
2598
2599 /* ------------------------------------------------------------------ */
2600 /* decFloatNextPlus -- next towards +Infinity                         */
2601 /*                                                                    */
2602 /*   result gets the next larger decFloat                             */
2603 /*   dfl    is the decFloat to start with                             */
2604 /*   set    is the context                                            */
2605 /*   returns result                                                   */
2606 /*                                                                    */
2607 /* This is 754r nextup; Invalid is the only status possible (from     */
2608 /* an sNaN).                                                          */
2609 /* ------------------------------------------------------------------ */
2610 decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
2611                             decContext *set) {
2612   uInt savestat;                        /* saves status */
2613   enum rounding saveround;              /* .. and mode */
2614   decFloat delta;                       /* tiny increment */
2615
2616   /* -Infinity is the special case */
2617   if (DFISINF(dfl) && DFISSIGNED(dfl)) {
2618     DFSETNMAX(result);
2619     DFWORD(result, 0)|=DECFLOAT_Sign;   /* make negative */
2620     return result;                      /* [no status to set] */
2621     }
2622   /* other cases are effected by sutracting a tiny delta -- this */
2623   /* should be done in a wider format as the delta is unrepresentable */
2624   /* here (but can be done with normal add if the sign of zero is */
2625   /* treated carefully, because no Inexactitude is interesting); */
2626   /* rounding to +Infinity then pushes the result to next above */
2627   decFloatZero(&delta);                 /* set up tiny delta */
2628   DFWORD(&delta, DECWORDS-1)=1;         /* coefficient=1 */
2629   DFWORD(&delta, 0)=0;                  /* Sign=0 + biased exponent=0 */
2630   /* set up for the directional round */
2631   saveround=set->round;                 /* save mode */
2632   set->round=DEC_ROUND_CEILING;         /* .. round towards +Infinity */
2633   savestat=set->status;                 /* save status */
2634   decFloatAdd(result, dfl, &delta, set);
2635   /* Add rules mess up the sign when going from -Ntiny to -0 */
2636   if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
2637   set->status&=DEC_Invalid_operation;   /* preserve only sNaN status */
2638   set->status|=savestat;                /* restore pending flags */
2639   set->round=saveround;                 /* .. and mode */
2640   return result;
2641   } /* decFloatNextPlus */
2642
2643 /* ------------------------------------------------------------------ */
2644 /* decFloatNextToward -- next towards a decFloat                      */
2645 /*                                                                    */
2646 /*   result gets the next decFloat                                    */
2647 /*   dfl    is the decFloat to start with                             */
2648 /*   dfr    is the decFloat to move toward                            */
2649 /*   set    is the context                                            */
2650 /*   returns result                                                   */
2651 /*                                                                    */
2652 /* This is 754r nextafter; status may be set unless the result is a   */
2653 /* normal number.                                                     */
2654 /* ------------------------------------------------------------------ */
2655 decFloat * decFloatNextToward(decFloat *result,
2656                               const decFloat *dfl, const decFloat *dfr,
2657                               decContext *set) {
2658   decFloat delta;                       /* tiny increment or decrement */
2659   decFloat pointone;                    /* 1e-1 */
2660   uInt  savestat;                       /* saves status */
2661   enum  rounding saveround;             /* .. and mode */
2662   uInt  deltatop;                       /* top word for delta */
2663   Int   comp;                           /* work */
2664
2665   if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2666   /* Both are numeric, so Invalid no longer a possibility */
2667   comp=decNumCompare(dfl, dfr, 0);
2668   if (comp==0) return decFloatCopySign(result, dfl, dfr); /* equal */
2669   /* unequal; do NextPlus or NextMinus but with different status rules */
2670
2671   if (comp<0) { /* lhs<rhs, do NextPlus, see above for commentary */
2672     if (DFISINF(dfl) && DFISSIGNED(dfl)) {   /* -Infinity special case */
2673       DFSETNMAX(result);
2674       DFWORD(result, 0)|=DECFLOAT_Sign;
2675       return result;
2676       }
2677     saveround=set->round;                    /* save mode */
2678     set->round=DEC_ROUND_CEILING;            /* .. round towards +Infinity */
2679     deltatop=0;                              /* positive delta */
2680     }
2681    else { /* lhs>rhs, do NextMinus, see above for commentary */
2682     if (DFISINF(dfl) && !DFISSIGNED(dfl)) {  /* +Infinity special case */
2683       DFSETNMAX(result);
2684       return result;
2685       }
2686     saveround=set->round;                    /* save mode */
2687     set->round=DEC_ROUND_FLOOR;              /* .. round towards -Infinity */
2688     deltatop=DECFLOAT_Sign;                  /* negative delta */
2689     }
2690   savestat=set->status;                      /* save status */
2691   /* Here, Inexact is needed where appropriate (and hence Underflow, */
2692   /* etc.).  Therefore the tiny delta which is otherwise */
2693   /* unrepresentable (see NextPlus and NextMinus) is constructed */
2694   /* using the multiplication of FMA. */
2695   decFloatZero(&delta);                 /* set up tiny delta */
2696   DFWORD(&delta, DECWORDS-1)=1;         /* coefficient=1 */
2697   DFWORD(&delta, 0)=deltatop;           /* Sign + biased exponent=0 */
2698   decFloatFromString(&pointone, "1E-1", set); /* set up multiplier */
2699   decFloatFMA(result, &delta, &pointone, dfl, set);
2700   /* [Delta is truly tiny, so no need to correct sign of zero] */
2701   /* use new status unless the result is normal */
2702   if (decFloatIsNormal(result)) set->status=savestat; /* else goes forward */
2703   set->round=saveround;                 /* restore mode */
2704   return result;
2705   } /* decFloatNextToward */
2706
2707 /* ------------------------------------------------------------------ */
2708 /* decFloatOr -- logical digitwise OR of two decFloats                */
2709 /*                                                                    */
2710 /*   result gets the result of ORing dfl and dfr                      */
2711 /*   dfl    is the first decFloat (lhs)                               */
2712 /*   dfr    is the second decFloat (rhs)                              */
2713 /*   set    is the context                                            */
2714 /*   returns result, which will be canonical with sign=0              */
2715 /*                                                                    */
2716 /* The operands must be positive, finite with exponent q=0, and       */
2717 /* comprise just zeros and ones; if not, Invalid operation results.   */
2718 /* ------------------------------------------------------------------ */
2719 decFloat * decFloatOr(decFloat *result,
2720                        const decFloat *dfl, const decFloat *dfr,
2721                        decContext *set) {
2722   if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
2723    || !DFISCC01(dfl)   || !DFISCC01(dfr)) return decInvalid(result, set);
2724   /* the operands are positive finite integers (q=0) with just 0s and 1s */
2725   #if DOUBLE
2726    DFWORD(result, 0)=ZEROWORD
2727                    |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04009124);
2728    DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x49124491;
2729   #elif QUAD
2730    DFWORD(result, 0)=ZEROWORD
2731                    |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04000912);
2732    DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x44912449;
2733    DFWORD(result, 2)=(DFWORD(dfl, 2) | DFWORD(dfr, 2))&0x12449124;
2734    DFWORD(result, 3)=(DFWORD(dfl, 3) | DFWORD(dfr, 3))&0x49124491;
2735   #endif
2736   return result;
2737   } /* decFloatOr */
2738
2739 /* ------------------------------------------------------------------ */
2740 /* decFloatPlus -- add value to 0, heeding NaNs, etc.                 */
2741 /*                                                                    */
2742 /*   result gets the canonicalized 0+df                               */
2743 /*   df     is the decFloat to plus                                   */
2744 /*   set    is the context                                            */
2745 /*   returns result                                                   */
2746 /*                                                                    */
2747 /* This has the same effect as 0+df where the exponent of the zero is */
2748 /* the same as that of df (if df is finite).                          */
2749 /* The effect is also the same as decFloatCopy except that NaNs       */
2750 /* are handled normally (the sign of a NaN is not affected, and an    */
2751 /* sNaN will signal), the result is canonical, and zero gets sign 0.  */
2752 /* ------------------------------------------------------------------ */
2753 decFloat * decFloatPlus(decFloat *result, const decFloat *df,
2754                         decContext *set) {
2755   if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2756   decCanonical(result, df);                       /* copy and check */
2757   if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80;     /* turn off sign bit */
2758   return result;
2759   } /* decFloatPlus */
2760
2761 /* ------------------------------------------------------------------ */
2762 /* decFloatQuantize -- quantize a decFloat                            */
2763 /*                                                                    */
2764 /*   result gets the result of quantizing dfl to match dfr            */
2765 /*   dfl    is the first decFloat (lhs)                               */
2766 /*   dfr    is the second decFloat (rhs), which sets the exponent     */
2767 /*   set    is the context                                            */
2768 /*   returns result                                                   */
2769 /*                                                                    */
2770 /* Unless there is an error or the result is infinite, the exponent   */
2771 /* of result is guaranteed to be the same as that of dfr.             */
2772 /* ------------------------------------------------------------------ */
2773 decFloat * decFloatQuantize(decFloat *result,
2774                             const decFloat *dfl, const decFloat *dfr,
2775                             decContext *set) {
2776   Int   explb, exprb;         /* left and right biased exponents */
2777   uByte *ulsd;                /* local LSD pointer */
2778   uInt  *ui;                  /* work */
2779   uByte *ub;                  /* .. */
2780   Int   drop;                 /* .. */
2781   uInt  dpd;                  /* .. */
2782   uInt  encode;               /* encoding accumulator */
2783   uInt  sourhil, sourhir;     /* top words from source decFloats */
2784   /* the following buffer holds the coefficient for manipulation */
2785   uByte buf[4+DECPMAX*3];     /* + space for zeros to left or right */
2786   #if DECTRACE
2787   bcdnum num;                 /* for trace displays */
2788   #endif
2789
2790   /* Start decoding the arguments */
2791   sourhil=DFWORD(dfl, 0);          /* LHS top word */
2792   explb=DECCOMBEXP[sourhil>>26];   /* get exponent high bits (in place) */
2793   sourhir=DFWORD(dfr, 0);          /* RHS top word */
2794   exprb=DECCOMBEXP[sourhir>>26];
2795
2796   if (EXPISSPECIAL(explb | exprb)) { /* either is special? */
2797     /* NaNs are handled as usual */
2798     if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2799     /* one infinity but not both is bad */
2800     if (DFISINF(dfl)!=DFISINF(dfr)) return decInvalid(result, set);
2801     /* both infinite; return canonical infinity with sign of LHS */
2802     return decInfinity(result, dfl);
2803     }
2804
2805   /* Here when both arguments are finite */
2806   /* complete extraction of the exponents [no need to unbias] */
2807   explb+=GETECON(dfl);             /* + continuation */
2808   exprb+=GETECON(dfr);             /* .. */
2809
2810   /* calculate the number of digits to drop from the coefficient */
2811   drop=exprb-explb;                /* 0 if nothing to do */
2812   if (drop==0) return decCanonical(result, dfl); /* return canonical */
2813
2814   /* the coefficient is needed; lay it out into buf, offset so zeros */
2815   /* can be added before or after as needed -- an extra heading is */
2816   /* added so can safely pad Quad DECPMAX-1 zeros to the left by */
2817   /* fours */
2818   #define BUFOFF (buf+4+DECPMAX)
2819   GETCOEFF(dfl, BUFOFF);           /* decode from decFloat */
2820   /* [now the msd is at BUFOFF and the lsd is at BUFOFF+DECPMAX-1] */
2821
2822   #if DECTRACE
2823   num.msd=BUFOFF;
2824   num.lsd=BUFOFF+DECPMAX-1;
2825   num.exponent=explb-DECBIAS;
2826   num.sign=sourhil & DECFLOAT_Sign;
2827   decShowNum(&num, "dfl");
2828   #endif
2829
2830   if (drop>0) {                         /* [most common case] */
2831     /* (this code is very similar to that in decFloatFinalize, but */
2832     /* has many differences so is duplicated here -- so any changes */
2833     /* may need to be made there, too) */
2834     uByte *roundat;                          /* -> re-round digit */
2835     uByte reround;                           /* reround value */
2836     /* printf("Rounding; drop=%ld\n", (LI)drop); */
2837
2838     /* there is at least one zero needed to the left, in all but one */
2839     /* exceptional (all-nines) case, so place four zeros now; this is */
2840     /* needed almost always and makes rounding all-nines by fours safe */
2841     UINTAT(BUFOFF-4)=0;
2842
2843     /* Three cases here: */
2844     /*   1. new LSD is in coefficient (almost always) */
2845     /*   2. new LSD is digit to left of coefficient (so MSD is */
2846     /*      round-for-reround digit) */
2847     /*   3. new LSD is to left of case 2 (whole coefficient is sticky) */
2848     /* Note that leading zeros can safely be treated as useful digits */
2849
2850     /* [duplicate check-stickies code to save a test] */
2851     /* [by-digit check for stickies as runs of zeros are rare] */
2852     if (drop<DECPMAX) {                      /* NB lengths not addresses */
2853       roundat=BUFOFF+DECPMAX-drop;
2854       reround=*roundat;
2855       for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
2856         if (*ub!=0) {                        /* non-zero to be discarded */
2857           reround=DECSTICKYTAB[reround];     /* apply sticky bit */
2858           break;                             /* [remainder don't-care] */
2859           }
2860         } /* check stickies */
2861       ulsd=roundat-1;                        /* set LSD */
2862       }
2863      else {                                  /* edge case */
2864       if (drop==DECPMAX) {
2865         roundat=BUFOFF;
2866         reround=*roundat;
2867         }
2868        else {
2869         roundat=BUFOFF-1;
2870         reround=0;
2871         }
2872       for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
2873         if (*ub!=0) {                        /* non-zero to be discarded */
2874           reround=DECSTICKYTAB[reround];     /* apply sticky bit */
2875           break;                             /* [remainder don't-care] */
2876           }
2877         } /* check stickies */
2878       *BUFOFF=0;                             /* make a coefficient of 0 */
2879       ulsd=BUFOFF;                           /* .. at the MSD place */
2880       }
2881
2882     if (reround!=0) {                        /* discarding non-zero */
2883       uInt bump=0;
2884       set->status|=DEC_Inexact;
2885
2886       /* next decide whether to increment the coefficient */
2887       if (set->round==DEC_ROUND_HALF_EVEN) { /* fastpath slowest case */
2888         if (reround>5) bump=1;               /* >0.5 goes up */
2889          else if (reround==5)                /* exactly 0.5000 .. */
2890           bump=*ulsd & 0x01;                 /* .. up iff [new] lsd is odd */
2891         } /* r-h-e */
2892        else switch (set->round) {
2893         case DEC_ROUND_DOWN: {
2894           /* no change */
2895           break;} /* r-d */
2896         case DEC_ROUND_HALF_DOWN: {
2897           if (reround>5) bump=1;
2898           break;} /* r-h-d */
2899         case DEC_ROUND_HALF_UP: {
2900           if (reround>=5) bump=1;
2901           break;} /* r-h-u */
2902         case DEC_ROUND_UP: {
2903           if (reround>0) bump=1;
2904           break;} /* r-u */
2905         case DEC_ROUND_CEILING: {
2906           /* same as _UP for positive numbers, and as _DOWN for negatives */
2907           if (!(sourhil&DECFLOAT_Sign) && reround>0) bump=1;
2908           break;} /* r-c */
2909         case DEC_ROUND_FLOOR: {
2910           /* same as _UP for negative numbers, and as _DOWN for positive */
2911           /* [negative reround cannot occur on 0] */
2912           if (sourhil&DECFLOAT_Sign && reround>0) bump=1;
2913           break;} /* r-f */
2914         case DEC_ROUND_05UP: {
2915           if (reround>0) { /* anything out there is 'sticky' */
2916             /* bump iff lsd=0 or 5; this cannot carry so it could be */
2917             /* effected immediately with no bump -- but the code */
2918             /* is clearer if this is done the same way as the others */
2919             if (*ulsd==0 || *ulsd==5) bump=1;
2920             }
2921           break;} /* r-r */
2922         default: {      /* e.g., DEC_ROUND_MAX */
2923           set->status|=DEC_Invalid_context;
2924           #if DECCHECK
2925           printf("Unknown rounding mode: %ld\n", (LI)set->round);
2926           #endif
2927           break;}
2928         } /* switch (not r-h-e) */
2929       /* printf("ReRound: %ld  bump: %ld\n", (LI)reround, (LI)bump); */
2930
2931       if (bump!=0) {                         /* need increment */
2932         /* increment the coefficient; this could give 1000... (after */
2933         /* the all nines case) */
2934         ub=ulsd;
2935         for (; UINTAT(ub-3)==0x09090909; ub-=4) UINTAT(ub-3)=0;
2936         /* now at most 3 digits left to non-9 (usually just the one) */
2937         for (; *ub==9; ub--) *ub=0;
2938         *ub+=1;
2939         /* [the all-nines case will have carried one digit to the */
2940         /* left of the original MSD -- just where it is needed] */
2941         } /* bump needed */
2942       } /* inexact rounding */
2943
2944     /* now clear zeros to the left so exactly DECPMAX digits will be */
2945     /* available in the coefficent -- the first word to the left was */
2946     /* cleared earlier for safe carry; now add any more needed */
2947     if (drop>4) {
2948       UINTAT(BUFOFF-8)=0;                    /* must be at least 5 */
2949       for (ui=&UINTAT(BUFOFF-12); ui>&UINTAT(ulsd-DECPMAX-3); ui--) *ui=0;
2950       }
2951     } /* need round (drop>0) */
2952
2953    else { /* drop<0; padding with -drop digits is needed */
2954     /* This is the case where an error can occur if the padded */
2955     /* coefficient will not fit; checking for this can be done in the */
2956     /* same loop as padding for zeros if the no-hope and zero cases */
2957     /* are checked first */
2958     if (-drop>DECPMAX-1) {                   /* cannot fit unless 0 */
2959       if (!ISCOEFFZERO(BUFOFF)) return decInvalid(result, set);
2960       /* a zero can have any exponent; just drop through and use it */
2961       ulsd=BUFOFF+DECPMAX-1;
2962       }
2963      else { /* padding will fit (but may still be too long) */
2964       /* final-word mask depends on endianess */
2965       #if DECLITEND
2966       static const uInt dmask[]={0, 0x000000ff, 0x0000ffff, 0x00ffffff};
2967       #else
2968       static const uInt dmask[]={0, 0xff000000, 0xffff0000, 0xffffff00};
2969       #endif
2970       for (ui=&UINTAT(BUFOFF+DECPMAX);; ui++) {
2971         *ui=0;
2972         if (UINTAT(&UBYTEAT(ui)-DECPMAX)!=0) { /* could be bad */
2973           /* if all four digits should be zero, definitely bad */
2974           if (ui<=&UINTAT(BUFOFF+DECPMAX+(-drop)-4))
2975             return decInvalid(result, set);
2976           /* must be a 1- to 3-digit sequence; check more carefully */
2977           if ((UINTAT(&UBYTEAT(ui)-DECPMAX)&dmask[(-drop)%4])!=0)
2978             return decInvalid(result, set);
2979           break;    /* no need for loop end test */
2980           }
2981         if (ui>=&UINTAT(BUFOFF+DECPMAX+(-drop)-4)) break; /* done */
2982         }
2983       ulsd=BUFOFF+DECPMAX+(-drop)-1;
2984       } /* pad and check leading zeros */
2985     } /* drop<0 */
2986
2987   #if DECTRACE
2988   num.msd=ulsd-DECPMAX+1;
2989   num.lsd=ulsd;
2990   num.exponent=explb-DECBIAS;
2991   num.sign=sourhil & DECFLOAT_Sign;
2992   decShowNum(&num, "res");
2993   #endif
2994
2995   /*------------------------------------------------------------------*/
2996   /* At this point the result is DECPMAX digits, ending at ulsd, so   */
2997   /* fits the encoding exactly; there is no possibility of error      */
2998   /*------------------------------------------------------------------*/
2999   encode=((exprb>>DECECONL)<<4) + *(ulsd-DECPMAX+1); /* make index */
3000   encode=DECCOMBFROM[encode];                /* indexed by (0-2)*16+msd */
3001   /* the exponent continuation can be extracted from the original RHS */
3002   encode|=sourhir & ECONMASK;
3003   encode|=sourhil&DECFLOAT_Sign;             /* add the sign from LHS */
3004
3005   /* finally encode the coefficient */
3006   /* private macro to encode a declet; this version can be used */
3007   /* because all coefficient digits exist */
3008   #define getDPD3q(dpd, n) ub=ulsd-(3*(n))-2;                   \
3009     dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];
3010
3011   #if DOUBLE
3012     getDPD3q(dpd, 4); encode|=dpd<<8;
3013     getDPD3q(dpd, 3); encode|=dpd>>2;
3014     DFWORD(result, 0)=encode;
3015     encode=dpd<<30;
3016     getDPD3q(dpd, 2); encode|=dpd<<20;
3017     getDPD3q(dpd, 1); encode|=dpd<<10;
3018     getDPD3q(dpd, 0); encode|=dpd;
3019     DFWORD(result, 1)=encode;
3020
3021   #elif QUAD
3022     getDPD3q(dpd,10); encode|=dpd<<4;
3023     getDPD3q(dpd, 9); encode|=dpd>>6;
3024     DFWORD(result, 0)=encode;
3025     encode=dpd<<26;
3026     getDPD3q(dpd, 8); encode|=dpd<<16;
3027     getDPD3q(dpd, 7); encode|=dpd<<6;
3028     getDPD3q(dpd, 6); encode|=dpd>>4;
3029     DFWORD(result, 1)=encode;
3030     encode=dpd<<28;
3031     getDPD3q(dpd, 5); encode|=dpd<<18;
3032     getDPD3q(dpd, 4); encode|=dpd<<8;
3033     getDPD3q(dpd, 3); encode|=dpd>>2;
3034     DFWORD(result, 2)=encode;
3035     encode=dpd<<30;
3036     getDPD3q(dpd, 2); encode|=dpd<<20;
3037     getDPD3q(dpd, 1); encode|=dpd<<10;
3038     getDPD3q(dpd, 0); encode|=dpd;
3039     DFWORD(result, 3)=encode;
3040   #endif
3041   return result;
3042   } /* decFloatQuantize */
3043
3044 /* ------------------------------------------------------------------ */
3045 /* decFloatReduce -- reduce finite coefficient to minimum length      */
3046 /*                                                                    */
3047 /*   result gets the reduced decFloat                                 */
3048 /*   df     is the source decFloat                                    */
3049 /*   set    is the context                                            */
3050 /*   returns result, which will be canonical                          */
3051 /*                                                                    */
3052 /* This removes all possible trailing zeros from the coefficient;     */
3053 /* some may remain when the number is very close to Nmax.             */
3054 /* Special values are unchanged and no status is set unless df=sNaN.  */
3055 /* Reduced zero has an exponent q=0.                                  */
3056 /* ------------------------------------------------------------------ */
3057 decFloat * decFloatReduce(decFloat *result, const decFloat *df,
3058                           decContext *set) {
3059   bcdnum num;                           /* work */
3060   uByte buf[DECPMAX], *ub;              /* coefficient and pointer */
3061   if (df!=result) *result=*df;          /* copy, if needed */
3062   if (DFISNAN(df)) return decNaNs(result, df, NULL, set);   /* sNaN */
3063   /* zeros and infinites propagate too */
3064   if (DFISINF(df)) return decInfinity(result, df);     /* canonical */
3065   if (DFISZERO(df)) {
3066     uInt sign=DFWORD(df, 0)&DECFLOAT_Sign;
3067     decFloatZero(result);
3068     DFWORD(result, 0)|=sign;
3069     return result;                      /* exponent dropped, sign OK */
3070     }
3071   /* non-zero finite */
3072   GETCOEFF(df, buf);
3073   ub=buf+DECPMAX-1;                     /* -> lsd */
3074   if (*ub) return result;               /* no trailing zeros */
3075   for (ub--; *ub==0;) ub--;             /* terminates because non-zero */
3076   /* *ub is the first non-zero from the right */
3077   num.sign=DFWORD(df, 0)&DECFLOAT_Sign; /* set up number... */
3078   num.exponent=GETEXPUN(df)+(Int)(buf+DECPMAX-1-ub); /* adjusted exponent */
3079   num.msd=buf;
3080   num.lsd=ub;
3081   return decFinalize(result, &num, set);
3082   } /* decFloatReduce */
3083
3084 /* ------------------------------------------------------------------ */
3085 /* decFloatRemainder -- integer divide and return remainder           */
3086 /*                                                                    */
3087 /*   result gets the remainder of dividing dfl by dfr:                */
3088 /*   dfl    is the first decFloat (lhs)                               */
3089 /*   dfr    is the second decFloat (rhs)                              */
3090 /*   set    is the context                                            */
3091 /*   returns result                                                   */
3092 /*                                                                    */
3093 /* ------------------------------------------------------------------ */
3094 decFloat * decFloatRemainder(decFloat *result,
3095                              const decFloat *dfl, const decFloat *dfr,
3096                              decContext *set) {
3097   return decDivide(result, dfl, dfr, set, REMAINDER);
3098   } /* decFloatRemainder */
3099
3100 /* ------------------------------------------------------------------ */
3101 /* decFloatRemainderNear -- integer divide to nearest and remainder   */
3102 /*                                                                    */
3103 /*   result gets the remainder of dividing dfl by dfr:                */
3104 /*   dfl    is the first decFloat (lhs)                               */
3105 /*   dfr    is the second decFloat (rhs)                              */
3106 /*   set    is the context                                            */
3107 /*   returns result                                                   */
3108 /*                                                                    */
3109 /* This is the IEEE remainder, where the nearest integer is used.     */
3110 /* ------------------------------------------------------------------ */
3111 decFloat * decFloatRemainderNear(decFloat *result,
3112                              const decFloat *dfl, const decFloat *dfr,
3113                              decContext *set) {
3114   return decDivide(result, dfl, dfr, set, REMNEAR);
3115   } /* decFloatRemainderNear */
3116
3117 /* ------------------------------------------------------------------ */
3118 /* decFloatRotate -- rotate the coefficient of a decFloat left/right  */
3119 /*                                                                    */
3120 /*   result gets the result of rotating dfl                           */
3121 /*   dfl    is the source decFloat to rotate                          */
3122 /*   dfr    is the count of digits to rotate, an integer (with q=0)   */
3123 /*   set    is the context                                            */
3124 /*   returns result                                                   */
3125 /*                                                                    */
3126 /* The digits of the coefficient of dfl are rotated to the left (if   */
3127 /* dfr is positive) or to the right (if dfr is negative) without      */
3128 /* adjusting the exponent or the sign of dfl.                         */
3129 /*                                                                    */
3130 /* dfr must be in the range -DECPMAX through +DECPMAX.                */
3131 /* NaNs are propagated as usual.  An infinite dfl is unaffected (but  */
3132 /* dfr must be valid).  No status is set unless dfr is invalid or an  */
3133 /* operand is an sNaN.  The result is canonical.                      */
3134 /* ------------------------------------------------------------------ */
3135 #define PHALF (ROUNDUP(DECPMAX/2, 4))   /* half length, rounded up */
3136 decFloat * decFloatRotate(decFloat *result,
3137                          const decFloat *dfl, const decFloat *dfr,
3138                          decContext *set) {
3139   Int rotate;                           /* dfr as an Int */
3140   uByte buf[DECPMAX+PHALF];             /* coefficient + half */
3141   uInt digits, savestat;                /* work */
3142   bcdnum num;                           /* .. */
3143   uByte *ub;                            /* .. */
3144
3145   if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3146   if (!DFISINT(dfr)) return decInvalid(result, set);
3147   digits=decFloatDigits(dfr);                    /* calculate digits */
3148   if (digits>2) return decInvalid(result, set);  /* definitely out of range */
3149   rotate=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
3150   if (rotate>DECPMAX) return decInvalid(result, set); /* too big */
3151   /* [from here on no error or status change is possible] */
3152   if (DFISINF(dfl)) return decInfinity(result, dfl);  /* canonical */
3153   /* handle no-rotate cases */
3154   if (rotate==0 || rotate==DECPMAX) return decCanonical(result, dfl);
3155   /* a real rotate is needed: 0 < rotate < DECPMAX */
3156   /* reduce the rotation to no more than half to reduce copying later */
3157   /* (for QUAD in fact half + 2 digits) */
3158   if (DFISSIGNED(dfr)) rotate=-rotate;
3159   if (abs(rotate)>PHALF) {
3160     if (rotate<0) rotate=DECPMAX+rotate;
3161      else rotate=rotate-DECPMAX;
3162     }
3163   /* now lay out the coefficient, leaving room to the right or the */
3164   /* left depending on the direction of rotation */
3165   ub=buf;
3166   if (rotate<0) ub+=PHALF;    /* rotate right, so space to left */
3167   GETCOEFF(dfl, ub);
3168   /* copy half the digits to left or right, and set num.msd */
3169   if (rotate<0) {
3170     memcpy(buf, buf+DECPMAX, PHALF);
3171     num.msd=buf+PHALF+rotate;
3172     }
3173    else {
3174     memcpy(buf+DECPMAX, buf, PHALF);
3175     num.msd=buf+rotate;
3176     }
3177   /* fill in rest of num */
3178   num.lsd=num.msd+DECPMAX-1;
3179   num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
3180   num.exponent=GETEXPUN(dfl);
3181   savestat=set->status;                 /* record */
3182   decFinalize(result, &num, set);
3183   set->status=savestat;                 /* restore */
3184   return result;
3185   } /* decFloatRotate */
3186
3187 /* ------------------------------------------------------------------ */
3188 /* decFloatSameQuantum -- test decFloats for same quantum             */
3189 /*                                                                    */
3190 /*   dfl    is the first decFloat (lhs)                               */
3191 /*   dfr    is the second decFloat (rhs)                              */
3192 /*   returns 1 if the operands have the same quantum, 0 otherwise     */
3193 /*                                                                    */
3194 /* No error is possible and no status results.                        */
3195 /* ------------------------------------------------------------------ */
3196 uInt decFloatSameQuantum(const decFloat *dfl, const decFloat *dfr) {
3197   if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) {
3198     if (DFISNAN(dfl) && DFISNAN(dfr)) return 1;
3199     if (DFISINF(dfl) && DFISINF(dfr)) return 1;
3200     return 0;  /* any other special mixture gives false */
3201     }
3202   if (GETEXP(dfl)==GETEXP(dfr)) return 1; /* biased exponents match */
3203   return 0;
3204   } /* decFloatSameQuantum */
3205
3206 /* ------------------------------------------------------------------ */
3207 /* decFloatScaleB -- multiply by a power of 10, as per 754r           */
3208 /*                                                                    */
3209 /*   result gets the result of the operation                          */
3210 /*   dfl    is the first decFloat (lhs)                               */
3211 /*   dfr    is the second decFloat (rhs), am integer (with q=0)       */
3212 /*   set    is the context                                            */
3213 /*   returns result                                                   */
3214 /*                                                                    */
3215 /* This computes result=dfl x 10**dfr where dfr is an integer in the  */
3216 /* range +/-2*(emax+pmax), typically resulting from LogB.             */
3217 /* Underflow and Overflow (with Inexact) may occur.  NaNs propagate   */
3218 /* as usual.                                                          */
3219 /* ------------------------------------------------------------------ */
3220 #define SCALEBMAX 2*(DECEMAX+DECPMAX)   /* D=800, Q=12356 */
3221 decFloat * decFloatScaleB(decFloat *result,
3222                           const decFloat *dfl, const decFloat *dfr,
3223                           decContext *set) {
3224   uInt digits;                          /* work */
3225   Int  expr;                            /* dfr as an Int */
3226
3227   if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3228   if (!DFISINT(dfr)) return decInvalid(result, set);
3229   digits=decFloatDigits(dfr);                /* calculate digits */
3230
3231   #if DOUBLE
3232   if (digits>3) return decInvalid(result, set);   /* definitely out of range */
3233   expr=DPD2BIN[DFWORD(dfr, 1)&0x3ff];             /* must be in bottom declet */
3234   #elif QUAD
3235   if (digits>5) return decInvalid(result, set);   /* definitely out of range */
3236   expr=DPD2BIN[DFWORD(dfr, 3)&0x3ff]              /* in bottom 2 declets .. */
3237       +DPD2BIN[(DFWORD(dfr, 3)>>10)&0x3ff]*1000;  /* .. */
3238   #endif
3239   if (expr>SCALEBMAX) return decInvalid(result, set);  /* oops */
3240   /* [from now on no error possible] */
3241   if (DFISINF(dfl)) return decInfinity(result, dfl);   /* canonical */
3242   if (DFISSIGNED(dfr)) expr=-expr;
3243   /* dfl is finite and expr is valid */
3244   *result=*dfl;                              /* copy to target */
3245   return decFloatSetExponent(result, set, GETEXPUN(result)+expr);
3246   } /* decFloatScaleB */
3247
3248 /* ------------------------------------------------------------------ */
3249 /* decFloatShift -- shift the coefficient of a decFloat left or right */
3250 /*                                                                    */
3251 /*   result gets the result of shifting dfl                           */
3252 /*   dfl    is the source decFloat to shift                           */
3253 /*   dfr    is the count of digits to shift, an integer (with q=0)    */
3254 /*   set    is the context                                            */
3255 /*   returns result                                                   */
3256 /*                                                                    */
3257 /* The digits of the coefficient of dfl are shifted to the left (if   */
3258 /* dfr is positive) or to the right (if dfr is negative) without      */
3259 /* adjusting the exponent or the sign of dfl.                         */
3260 /*                                                                    */
3261 /* dfr must be in the range -DECPMAX through +DECPMAX.                */
3262 /* NaNs are propagated as usual.  An infinite dfl is unaffected (but  */
3263 /* dfr must be valid).  No status is set unless dfr is invalid or an  */
3264 /* operand is an sNaN.  The result is canonical.                      */
3265 /* ------------------------------------------------------------------ */
3266 decFloat * decFloatShift(decFloat *result,
3267                          const decFloat *dfl, const decFloat *dfr,
3268                          decContext *set) {
3269   Int shift;                            /* dfr as an Int */
3270   uByte buf[DECPMAX*2];                 /* coefficient + padding */
3271   uInt digits, savestat;                /* work */
3272   bcdnum num;                           /* .. */
3273
3274   if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3275   if (!DFISINT(dfr)) return decInvalid(result, set);
3276   digits=decFloatDigits(dfr);                     /* calculate digits */
3277   if (digits>2) return decInvalid(result, set);   /* definitely out of range */
3278   shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff];   /* is in bottom declet */
3279   if (shift>DECPMAX) return decInvalid(result, set);   /* too big */
3280   /* [from here on no error or status change is possible] */
3281
3282   if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
3283   /* handle no-shift and all-shift (clear to zero) cases */
3284   if (shift==0) return decCanonical(result, dfl);
3285   if (shift==DECPMAX) {                      /* zero with sign */
3286     uByte sign=(uByte)(DFBYTE(dfl, 0)&0x80); /* save sign bit */
3287     decFloatZero(result);                    /* make +0 */
3288     DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* and set sign */
3289     /* [cannot safely use CopySign] */
3290     return result;
3291     }
3292   /* a real shift is needed: 0 < shift < DECPMAX */
3293   num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
3294   num.exponent=GETEXPUN(dfl);
3295   num.msd=buf;
3296   GETCOEFF(dfl, buf);
3297   if (DFISSIGNED(dfr)) { /* shift right */
3298     /* edge cases are taken care of, so this is easy */
3299     num.lsd=buf+DECPMAX-shift-1;
3300     }
3301    else { /* shift left -- zero padding needed to right */
3302     UINTAT(buf+DECPMAX)=0;              /* 8 will handle most cases */
3303     UINTAT(buf+DECPMAX+4)=0;            /* .. */
3304     if (shift>8) memset(buf+DECPMAX+8, 0, 8+QUAD*18); /* all other cases */
3305     num.msd+=shift;
3306     num.lsd=num.msd+DECPMAX-1;
3307     }
3308   savestat=set->status;                 /* record */
3309   decFinalize(result, &num, set);
3310   set->status=savestat;                 /* restore */
3311   return result;
3312   } /* decFloatShift */
3313
3314 /* ------------------------------------------------------------------ */
3315 /* decFloatSubtract -- subtract a decFloat from another               */
3316 /*                                                                    */
3317 /*   result gets the result of subtracting dfr from dfl:              */
3318 /*   dfl    is the first decFloat (lhs)                               */
3319 /*   dfr    is the second decFloat (rhs)                              */
3320 /*   set    is the context                                            */
3321 /*   returns result                                                   */
3322 /*                                                                    */
3323 /* ------------------------------------------------------------------ */
3324 decFloat * decFloatSubtract(decFloat *result,
3325                             const decFloat *dfl, const decFloat *dfr,
3326                             decContext *set) {
3327   decFloat temp;
3328   /* NaNs must propagate without sign change */
3329   if (DFISNAN(dfr)) return decFloatAdd(result, dfl, dfr, set);
3330   temp=*dfr;                                   /* make a copy */
3331   DFBYTE(&temp, 0)^=0x80;                      /* flip sign */
3332   return decFloatAdd(result, dfl, &temp, set); /* and add to the lhs */
3333   } /* decFloatSubtract */
3334
3335 /* ------------------------------------------------------------------ */
3336 /* decFloatToInt -- round to 32-bit binary integer (4 flavours)       */
3337 /*                                                                    */
3338 /*   df    is the decFloat to round                                   */
3339 /*   set   is the context                                             */
3340 /*   round is the rounding mode to use                                */
3341 /*   returns a uInt or an Int, rounded according to the name          */
3342 /*                                                                    */
3343 /* Invalid will always be signaled if df is a NaN, is Infinite, or is */
3344 /* outside the range of the target; Inexact will not be signaled for  */
3345 /* simple rounding unless 'Exact' appears in the name.                */
3346 /* ------------------------------------------------------------------ */
3347 uInt decFloatToUInt32(const decFloat *df, decContext *set,
3348                       enum rounding round) {
3349   return decToInt32(df, set, round, 0, 1);}
3350
3351 uInt decFloatToUInt32Exact(const decFloat *df, decContext *set,
3352                            enum rounding round) {
3353   return decToInt32(df, set, round, 1, 1);}
3354
3355 Int decFloatToInt32(const decFloat *df, decContext *set,
3356                     enum rounding round) {
3357   return (Int)decToInt32(df, set, round, 0, 0);}
3358
3359 Int decFloatToInt32Exact(const decFloat *df, decContext *set,
3360                          enum rounding round) {
3361   return (Int)decToInt32(df, set, round, 1, 0);}
3362
3363 /* ------------------------------------------------------------------ */
3364 /* decFloatToIntegral -- round to integral value (two flavours)       */
3365 /*                                                                    */
3366 /*   result gets the result                                           */
3367 /*   df     is the decFloat to round                                  */
3368 /*   set    is the context                                            */
3369 /*   round  is the rounding mode to use                               */
3370 /*   returns result                                                   */
3371 /*                                                                    */
3372 /* No exceptions, even Inexact, are raised except for sNaN input, or  */
3373 /* if 'Exact' appears in the name.                                    */
3374 /* ------------------------------------------------------------------ */
3375 decFloat * decFloatToIntegralValue(decFloat *result, const decFloat *df,
3376                                    decContext *set, enum rounding round) {
3377   return decToIntegral(result, df, set, round, 0);}
3378
3379 decFloat * decFloatToIntegralExact(decFloat *result, const decFloat *df,
3380                                    decContext *set) {
3381   return decToIntegral(result, df, set, set->round, 1);}
3382
3383 /* ------------------------------------------------------------------ */
3384 /* decFloatXor -- logical digitwise XOR of two decFloats              */
3385 /*                                                                    */
3386 /*   result gets the result of XORing dfl and dfr                     */
3387 /*   dfl    is the first decFloat (lhs)                               */
3388 /*   dfr    is the second decFloat (rhs)                              */
3389 /*   set    is the context                                            */
3390 /*   returns result, which will be canonical with sign=0              */
3391 /*                                                                    */
3392 /* The operands must be positive, finite with exponent q=0, and       */
3393 /* comprise just zeros and ones; if not, Invalid operation results.   */
3394 /* ------------------------------------------------------------------ */
3395 decFloat * decFloatXor(decFloat *result,
3396                        const decFloat *dfl, const decFloat *dfr,
3397                        decContext *set) {
3398   if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
3399    || !DFISCC01(dfl)   || !DFISCC01(dfr)) return decInvalid(result, set);
3400   /* the operands are positive finite integers (q=0) with just 0s and 1s */
3401   #if DOUBLE
3402    DFWORD(result, 0)=ZEROWORD
3403                    |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04009124);
3404    DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x49124491;
3405   #elif QUAD
3406    DFWORD(result, 0)=ZEROWORD
3407                    |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04000912);
3408    DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x44912449;
3409    DFWORD(result, 2)=(DFWORD(dfl, 2) ^ DFWORD(dfr, 2))&0x12449124;
3410    DFWORD(result, 3)=(DFWORD(dfl, 3) ^ DFWORD(dfr, 3))&0x49124491;
3411   #endif
3412   return result;
3413   } /* decFloatXor */
3414
3415 /* ------------------------------------------------------------------ */
3416 /* decInvalid -- set Invalid_operation result                         */
3417 /*                                                                    */
3418 /*   result gets a canonical NaN                                      */
3419 /*   set    is the context                                            */
3420 /*   returns result                                                   */
3421 /*                                                                    */
3422 /* status has Invalid_operation added                                 */
3423 /* ------------------------------------------------------------------ */
3424 static decFloat *decInvalid(decFloat *result, decContext *set) {
3425   decFloatZero(result);
3426   DFWORD(result, 0)=DECFLOAT_qNaN;
3427   set->status|=DEC_Invalid_operation;
3428   return result;
3429   } /* decInvalid */
3430
3431 /* ------------------------------------------------------------------ */
3432 /* decInfinity -- set canonical Infinity with sign from a decFloat    */
3433 /*                                                                    */
3434 /*   result gets a canonical Infinity                                 */
3435 /*   df     is source decFloat (only the sign is used)                */
3436 /*   returns result                                                   */
3437 /*                                                                    */
3438 /* df may be the same as result                                       */
3439 /* ------------------------------------------------------------------ */
3440 static decFloat *decInfinity(decFloat *result, const decFloat *df) {
3441   uInt sign=DFWORD(df, 0);         /* save source signword */
3442   decFloatZero(result);            /* clear everything */
3443   DFWORD(result, 0)=DECFLOAT_Inf | (sign & DECFLOAT_Sign);
3444   return result;
3445   } /* decInfinity */
3446
3447 /* ------------------------------------------------------------------ */
3448 /* decNaNs -- handle NaN argument(s)                                  */
3449 /*                                                                    */
3450 /*   result gets the result of handling dfl and dfr, one or both of   */
3451 /*          which is a NaN                                            */
3452 /*   dfl    is the first decFloat (lhs)                               */
3453 /*   dfr    is the second decFloat (rhs) -- may be NULL for a single- */
3454 /*          operand operation                                         */
3455 /*   set    is the context                                            */
3456 /*   returns result                                                   */
3457 /*                                                                    */
3458 /* Called when one or both operands is a NaN, and propagates the      */
3459 /* appropriate result to res.  When an sNaN is found, it is changed   */
3460 /* to a qNaN and Invalid operation is set.                            */
3461 /* ------------------------------------------------------------------ */
3462 static decFloat *decNaNs(decFloat *result,
3463                          const decFloat *dfl, const decFloat *dfr,
3464                          decContext *set) {
3465   /* handle sNaNs first */
3466   if (dfr!=NULL && DFISSNAN(dfr) && !DFISSNAN(dfl)) dfl=dfr; /* use RHS */
3467   if (DFISSNAN(dfl)) {
3468     decCanonical(result, dfl);          /* propagate canonical sNaN */
3469     DFWORD(result, 0)&=~(DECFLOAT_qNaN ^ DECFLOAT_sNaN); /* quiet */
3470     set->status|=DEC_Invalid_operation;
3471     return result;
3472     }
3473   /* one or both is a quiet NaN */
3474   if (!DFISNAN(dfl)) dfl=dfr;           /* RHS must be NaN, use it */
3475   return decCanonical(result, dfl);     /* propagate canonical qNaN */
3476   } /* decNaNs */
3477
3478 /* ------------------------------------------------------------------ */
3479 /* decNumCompare -- numeric comparison of two decFloats               */
3480 /*                                                                    */
3481 /*   dfl    is the left-hand decFloat, which is not a NaN             */
3482 /*   dfr    is the right-hand decFloat, which is not a NaN            */
3483 /*   tot    is 1 for total order compare, 0 for simple numeric        */
3484 /*   returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr               */
3485 /*                                                                    */
3486 /* No error is possible; status and mode are unchanged.               */
3487 /* ------------------------------------------------------------------ */
3488 static Int decNumCompare(const decFloat *dfl, const decFloat *dfr, Flag tot) {
3489   Int   sigl, sigr;                     /* LHS and RHS non-0 signums */
3490   Int   shift;                          /* shift needed to align operands */
3491   uByte *ub, *uc;                       /* work */
3492   /* buffers +2 if Quad (36 digits), need double plus 4 for safe padding */
3493   uByte bufl[DECPMAX*2+QUAD*2+4];       /* for LHS coefficient + padding */
3494   uByte bufr[DECPMAX*2+QUAD*2+4];       /* for RHS coefficient + padding */
3495
3496   sigl=1;
3497   if (DFISSIGNED(dfl)) {
3498     if (!DFISSIGNED(dfr)) {             /* -LHS +RHS */
3499       if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
3500       return -1;                        /* RHS wins */
3501       }
3502     sigl=-1;
3503     }
3504   if (DFISSIGNED(dfr)) {
3505     if (!DFISSIGNED(dfl)) {             /* +LHS -RHS */
3506       if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
3507       return +1;                        /* LHS wins */
3508       }
3509     }
3510
3511   /* signs are the same; operand(s) could be zero */
3512   sigr=-sigl;                           /* sign to return if abs(RHS) wins */
3513
3514   if (DFISINF(dfl)) {
3515     if (DFISINF(dfr)) return 0;         /* both infinite & same sign */
3516     return sigl;                        /* inf > n */
3517     }
3518   if (DFISINF(dfr)) return sigr;        /* n < inf [dfl is finite] */
3519
3520   /* here, both are same sign and finite; calculate their offset */
3521   shift=GETEXP(dfl)-GETEXP(dfr);        /* [0 means aligned] */
3522   /* [bias can be ignored -- the absolute exponent is not relevant] */
3523
3524   if (DFISZERO(dfl)) {
3525     if (!DFISZERO(dfr)) return sigr;    /* LHS=0, RHS!=0 */
3526     /* both are zero, return 0 if both same exponent or numeric compare */
3527     if (shift==0 || !tot) return 0;
3528     if (shift>0) return sigl;
3529     return sigr;                        /* [shift<0] */
3530     }
3531    else {                               /* LHS!=0 */
3532     if (DFISZERO(dfr)) return sigl;     /* LHS!=0, RHS=0 */
3533     }
3534   /* both are known to be non-zero at this point */
3535
3536   /* if the exponents are so different that the coefficients do not */
3537   /* overlap (by even one digit) then a full comparison is not needed */
3538   if (abs(shift)>=DECPMAX) {            /* no overlap */
3539     /* coefficients are known to be non-zero */
3540     if (shift>0) return sigl;
3541     return sigr;                        /* [shift<0] */
3542     }
3543
3544   /* decode the coefficients */
3545   /* (shift both right two if Quad to make a multiple of four) */
3546   #if QUAD
3547     ub=bufl;                            /* avoid type-pun violation */
3548     UINTAT(ub)=0;
3549     uc=bufr;                            /* avoid type-pun violation */
3550     UINTAT(uc)=0;
3551   #endif
3552   GETCOEFF(dfl, bufl+QUAD*2);           /* decode from decFloat */
3553   GETCOEFF(dfr, bufr+QUAD*2);           /* .. */
3554   if (shift==0) {                       /* aligned; common and easy */
3555     /* all multiples of four, here */
3556     for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
3557       if (UINTAT(ub)==UINTAT(uc)) continue; /* so far so same */
3558       /* about to find a winner; go by bytes in case little-endian */
3559       for (;; ub++, uc++) {
3560         if (*ub>*uc) return sigl;       /* difference found */
3561         if (*ub<*uc) return sigr;       /* .. */
3562         }
3563       }
3564     } /* aligned */
3565    else if (shift>0) {                  /* lhs to left */
3566     ub=bufl;                            /* RHS pointer */
3567     /* pad bufl so right-aligned; most shifts will fit in 8 */
3568     UINTAT(bufl+DECPMAX+QUAD*2)=0;      /* add eight zeros */
3569     UINTAT(bufl+DECPMAX+QUAD*2+4)=0;    /* .. */
3570     if (shift>8) {
3571       /* more than eight; fill the rest, and also worth doing the */
3572       /* lead-in by fours */
3573       uByte *up;                         /* work */
3574       uByte *upend=bufl+DECPMAX+QUAD*2+shift;
3575       for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
3576       /* [pads up to 36 in all for Quad] */
3577       for (;; ub+=4) {
3578         if (UINTAT(ub)!=0) return sigl;
3579         if (ub+4>bufl+shift-4) break;
3580         }
3581       }
3582     /* check remaining leading digits */
3583     for (; ub<bufl+shift; ub++) if (*ub!=0) return sigl;
3584     /* now start the overlapped part; bufl has been padded, so the */
3585     /* comparison can go for the full length of bufr, which is a */
3586     /* multiple of 4 bytes */
3587     for (uc=bufr; ; uc+=4, ub+=4) {
3588       if (UINTAT(uc)!=UINTAT(ub)) {     /* mismatch found */
3589         for (;; uc++, ub++) {           /* check from left [little-endian?] */
3590           if (*ub>*uc) return sigl;     /* difference found */
3591           if (*ub<*uc) return sigr;     /* .. */
3592           }
3593         } /* mismatch */
3594       if (uc==bufr+QUAD*2+DECPMAX-4) break; /* all checked */
3595       }
3596     } /* shift>0 */
3597
3598    else { /* shift<0) .. RHS is to left of LHS; mirror shift>0 */
3599     uc=bufr;                            /* RHS pointer */
3600     /* pad bufr so right-aligned; most shifts will fit in 8 */
3601     UINTAT(bufr+DECPMAX+QUAD*2)=0;      /* add eight zeros */
3602     UINTAT(bufr+DECPMAX+QUAD*2+4)=0;    /* .. */
3603     if (shift<-8) {
3604       /* more than eight; fill the rest, and also worth doing the */
3605       /* lead-in by fours */
3606       uByte *up;                         /* work */
3607       uByte *upend=bufr+DECPMAX+QUAD*2-shift;
3608       for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
3609       /* [pads up to 36 in all for Quad] */
3610       for (;; uc+=4) {
3611         if (UINTAT(uc)!=0) return sigr;
3612         if (uc+4>bufr-shift-4) break;
3613         }
3614       }
3615     /* check remaining leading digits */
3616     for (; uc<bufr-shift; uc++) if (*uc!=0) return sigr;
3617     /* now start the overlapped part; bufr has been padded, so the */
3618     /* comparison can go for the full length of bufl, which is a */
3619     /* multiple of 4 bytes */
3620     for (ub=bufl; ; ub+=4, uc+=4) {
3621       if (UINTAT(ub)!=UINTAT(uc)) {     /* mismatch found */
3622         for (;; ub++, uc++) {           /* check from left [little-endian?] */
3623           if (*ub>*uc) return sigl;     /* difference found */
3624           if (*ub<*uc) return sigr;     /* .. */
3625           }
3626         } /* mismatch */
3627       if (ub==bufl+QUAD*2+DECPMAX-4) break; /* all checked */
3628       }
3629     } /* shift<0 */
3630
3631   /* Here when compare equal */
3632   if (!tot) return 0;                   /* numerically equal */
3633   /* total ordering .. exponent matters */
3634   if (shift>0) return sigl;             /* total order by exponent */
3635   if (shift<0) return sigr;             /* .. */
3636   return 0;
3637   } /* decNumCompare */
3638
3639 /* ------------------------------------------------------------------ */
3640 /* decToInt32 -- local routine to effect ToInteger conversions        */
3641 /*                                                                    */
3642 /*   df     is the decFloat to convert                                */
3643 /*   set    is the context                                            */
3644 /*   rmode  is the rounding mode to use                               */
3645 /*   exact  is 1 if Inexact should be signalled                       */
3646 /*   unsign is 1 if the result a uInt, 0 if an Int (cast to uInt)     */
3647 /*   returns 32-bit result as a uInt                                  */
3648 /*                                                                    */
3649 /* Invalid is set is df is a NaN, is infinite, or is out-of-range; in */
3650 /* these cases 0 is returned.                                         */
3651 /* ------------------------------------------------------------------ */
3652 static uInt decToInt32(const decFloat *df, decContext *set,
3653                        enum rounding rmode, Flag exact, Flag unsign) {
3654   Int  exp;                        /* exponent */
3655   uInt sourhi, sourpen, sourlo;    /* top word from source decFloat .. */
3656   uInt hi, lo;                     /* .. penultimate, least, etc. */
3657   decFloat zero, result;           /* work */
3658   Int  i;                          /* .. */
3659
3660   /* Start decoding the argument */
3661   sourhi=DFWORD(df, 0);                 /* top word */
3662   exp=DECCOMBEXP[sourhi>>26];           /* get exponent high bits (in place) */
3663   if (EXPISSPECIAL(exp)) {              /* is special? */
3664     set->status|=DEC_Invalid_operation; /* signal */
3665     return 0;
3666     }
3667
3668   /* Here when the argument is finite */
3669   if (GETEXPUN(df)==0) result=*df;      /* already a true integer */
3670    else {                               /* need to round to integer */
3671     enum rounding saveround;            /* saver */
3672     uInt savestatus;                    /* .. */
3673     saveround=set->round;               /* save rounding mode .. */
3674     savestatus=set->status;             /* .. and status */
3675     set->round=rmode;                   /* set mode */
3676     decFloatZero(&zero);                /* make 0E+0 */
3677     set->status=0;                      /* clear */
3678     decFloatQuantize(&result, df, &zero, set); /* [this may fail] */
3679     set->round=saveround;               /* restore rounding mode .. */
3680     if (exact) set->status|=savestatus; /* include Inexact */
3681      else set->status=savestatus;       /* .. or just original status */
3682     }
3683
3684   /* only the last four declets of the coefficient can contain */
3685   /* non-zero; check for others (and also NaN or Infinity from the */
3686   /* Quantize) first (see DFISZERO for explanation): */
3687   /* decFloatShow(&result, "sofar"); */
3688   #if DOUBLE
3689   if ((DFWORD(&result, 0)&0x1c03ff00)!=0
3690    || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
3691   #elif QUAD
3692   if ((DFWORD(&result, 2)&0xffffff00)!=0
3693    ||  DFWORD(&result, 1)!=0
3694    || (DFWORD(&result, 0)&0x1c003fff)!=0
3695    || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
3696   #endif
3697     set->status|=DEC_Invalid_operation; /* Invalid or out of range */
3698     return 0;
3699     }
3700   /* get last twelve digits of the coefficent into hi & ho, base */
3701   /* 10**9 (see GETCOEFFBILL): */
3702   sourlo=DFWORD(&result, DECWORDS-1);
3703   lo=DPD2BIN0[sourlo&0x3ff]
3704     +DPD2BINK[(sourlo>>10)&0x3ff]
3705     +DPD2BINM[(sourlo>>20)&0x3ff];
3706   sourpen=DFWORD(&result, DECWORDS-2);
3707   hi=DPD2BIN0[((sourpen<<2) | (sourlo>>30))&0x3ff];
3708
3709   /* according to request, check range carefully */
3710   if (unsign) {
3711     if (hi>4 || (hi==4 && lo>294967295) || (hi+lo!=0 && DFISSIGNED(&result))) {
3712       set->status|=DEC_Invalid_operation; /* out of range */
3713       return 0;
3714       }
3715     return hi*BILLION+lo;
3716     }
3717   /* signed */
3718   if (hi>2 || (hi==2 && lo>147483647)) {
3719     /* handle the usual edge case */
3720     if (lo==147483648 && hi==2 && DFISSIGNED(&result)) return 0x80000000;
3721     set->status|=DEC_Invalid_operation; /* truly out of range */
3722     return 0;
3723     }
3724   i=hi*BILLION+lo;
3725   if (DFISSIGNED(&result)) i=-i;
3726   return (uInt)i;
3727   } /* decToInt32 */
3728
3729 /* ------------------------------------------------------------------ */
3730 /* decToIntegral -- local routine to effect ToIntegral value          */
3731 /*                                                                    */
3732 /*   result gets the result                                           */
3733 /*   df     is the decFloat to round                                  */
3734 /*   set    is the context                                            */
3735 /*   rmode  is the rounding mode to use                               */
3736 /*   exact  is 1 if Inexact should be signalled                       */
3737 /*   returns result                                                   */
3738 /* ------------------------------------------------------------------ */
3739 static decFloat * decToIntegral(decFloat *result, const decFloat *df,
3740                                 decContext *set, enum rounding rmode,
3741                                 Flag exact) {
3742   Int  exp;                        /* exponent */
3743   uInt sourhi;                     /* top word from source decFloat */
3744   enum rounding saveround;         /* saver */
3745   uInt savestatus;                 /* .. */
3746   decFloat zero;                   /* work */
3747
3748   /* Start decoding the argument */
3749   sourhi=DFWORD(df, 0);            /* top word */
3750   exp=DECCOMBEXP[sourhi>>26];      /* get exponent high bits (in place) */
3751
3752   if (EXPISSPECIAL(exp)) {         /* is special? */
3753     /* NaNs are handled as usual */
3754     if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
3755     /* must be infinite; return canonical infinity with sign of df */
3756     return decInfinity(result, df);
3757     }
3758
3759   /* Here when the argument is finite */
3760   /* complete extraction of the exponent */
3761   exp+=GETECON(df)-DECBIAS;             /* .. + continuation and unbias */
3762
3763   if (exp>=0) return decCanonical(result, df); /* already integral */
3764
3765   saveround=set->round;                 /* save rounding mode .. */
3766   savestatus=set->status;               /* .. and status */
3767   set->round=rmode;                     /* set mode */
3768   decFloatZero(&zero);                  /* make 0E+0 */
3769   decFloatQuantize(result, df, &zero, set); /* 'integrate'; cannot fail */
3770   set->round=saveround;                 /* restore rounding mode .. */
3771   if (!exact) set->status=savestatus;   /* .. and status, unless exact */
3772   return result;
3773   } /* decToIntegral */