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