text: /* Common base code for the decNumber C Library.
     Copyright (C) 2007, 2009 Free Software Foundation, Inc.
     Contributed by IBM Corporation.  Author Mike Cowlishaw.
  
     This file is part of GCC.
  
     GCC is free software; you can redistribute it and/or modify it under
     the terms of the GNU General Public License as published by the Free
     Software Foundation; either version 3, or (at your option) any later
     version.
  
     GCC is distributed in the hope that it will be useful, but WITHOUT ANY
     WARRANTY; without even the implied warranty of MERCHANTABILITY or
     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
     for more details.
  
  Under Section 7 of GPL version 3, you are granted additional
  permissions described in the GCC Runtime Library Exception, version
  3.1, as published by the Free Software Foundation.
  
  You should have received a copy of the GNU General Public License and
  a copy of the GCC Runtime Library Exception along with this program;
  see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  <http://www.gnu.org/licenses/>.  */
  
  /* ------------------------------------------------------------------ */
  /* decBasic.c -- common base code for Basic decimal types	      */
  /* ------------------------------------------------------------------ */
  /* This module comprises code that is shared between decDouble and    */
  /* decQuad (but not decSingle).	 The main arithmetic operations are   */
  /* here (Add, Subtract, Multiply, FMA, and Division operators).	      */
  /*								      */
  /* Unlike decNumber, parameterization takes place at compile time     */
  /* rather than at runtime.  The parameters are set in the decDouble.c */
  /* (etc.) files, which then include this one to produce the compiled  */
  /* code.  The functions here, therefore, are code shared between      */
  /* multiple formats.						      */
  /*								      */
  /* This must be included after decCommon.c.			      */
  /* ------------------------------------------------------------------ */
  /* Names here refer to decFloat rather than to decDouble, etc., and */
  /* the functions are in strict alphabetical order. */
  
  /* The compile-time flags SINGLE, DOUBLE, and QUAD are set up in */
  /* decCommon.c */
  #if !defined(QUAD)
    #error decBasic.c must be included after decCommon.c
  #endif
  #if SINGLE
    #error Routines in decBasic.c are for decDouble and decQuad only
  #endif
  
  /* Private constants */
  #define DIVIDE	    0x80000000	   /* Divide operations [as flags] */
  #define REMAINDER   0x40000000	   /* .. */
  #define DIVIDEINT   0x20000000	   /* .. */
  #define REMNEAR	    0x10000000	   /* .. */
  
  /* Private functions (local, used only by routines in this module) */
  static decFloat *decDivide(decFloat *, const decFloat *,
  			      const decFloat *, decContext *, uInt);
  static decFloat *decCanonical(decFloat *, const decFloat *);
  static void	 decFiniteMultiply(bcdnum *, uByte *, const decFloat *,
  			      const decFloat *);
  static decFloat *decInfinity(decFloat *, const decFloat *);
  static decFloat *decInvalid(decFloat *, decContext *);
  static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *,
  			      decContext *);
  static Int	 decNumCompare(const decFloat *, const decFloat *, Flag);
  static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *,
  			      enum rounding, Flag);
  static uInt	 decToInt32(const decFloat *, decContext *, enum rounding,
  			      Flag, Flag);
  
  /* ------------------------------------------------------------------ */
  /* decCanonical -- copy a decFloat, making canonical		      */
  /*								      */
  /*   result gets the canonicalized df				      */
  /*   df	    is the decFloat to copy and make canonical		      */
  /*   returns result						      */
  /*								      */
  /* This is exposed via decFloatCanonical for Double and Quad only.    */
  /* This works on specials, too; no error or exception is possible.    */
  /* ------------------------------------------------------------------ */
  static decFloat * decCanonical(decFloat *result, const decFloat *df) {
    uInt encode, precode, dpd;	   /* work */
    uInt inword, uoff, canon;	   /* .. */
    Int  n;			   /* counter (down) */
    if (df!=result) *result=*df;	   /* effect copy if needed */
    if (DFISSPECIAL(result)) {
      if (DFISINF(result)) return decInfinity(result, df); /* clean Infinity */
      /* is a NaN */
      DFWORD(result, 0)&=~ECONNANMASK;	/* clear ECON except selector */
      if (DFISCCZERO(df)) return result;	/* coefficient continuation is 0 */
      /* drop through to check payload */
      }
    /* return quickly if the coefficient continuation is canonical */
    { /* declare block */
    #if DOUBLE
      uInt sourhi=DFWORD(df, 0);
      uInt sourlo=DFWORD(df, 1);
      if (CANONDPDOFF(sourhi, 8)
       && CANONDPDTWO(sourhi, sourlo, 30)
       && CANONDPDOFF(sourlo, 20)
       && CANONDPDOFF(sourlo, 10)
       && CANONDPDOFF(sourlo, 0)) return result;
    #elif QUAD
      uInt sourhi=DFWORD(df, 0);
      uInt sourmh=DFWORD(df, 1);
      uInt sourml=DFWORD(df, 2);
      uInt sourlo=DFWORD(df, 3);
      if (CANONDPDOFF(sourhi, 4)
       && CANONDPDTWO(sourhi, sourmh, 26)
       && CANONDPDOFF(sourmh, 16)
       && CANONDPDOFF(sourmh, 6)
       && CANONDPDTWO(sourmh, sourml, 28)
       && CANONDPDOFF(sourml, 18)
       && CANONDPDOFF(sourml, 8)
       && CANONDPDTWO(sourml, sourlo, 30)
       && CANONDPDOFF(sourlo, 20)
       && CANONDPDOFF(sourlo, 10)
       && CANONDPDOFF(sourlo, 0)) return result;
    #endif
    } /* block */
  
    /* Loop to repair a non-canonical coefficent, as needed */
    inword=DECWORDS-1;		   /* current input word */
    uoff=0;			   /* bit offset of declet */
    encode=DFWORD(result, inword);
    for (n=DECLETS-1; n>=0; n--) {   /* count down declets of 10 bits */
      dpd=encode>>uoff;
      uoff+=10;
      if (uoff>32) {		   /* crossed uInt boundary */
        inword--;
        encode=DFWORD(result, inword);
        uoff-=32;
        dpd|=encode<<(10-uoff);	   /* get pending bits */
        }
      dpd&=0x3ff;			   /* clear uninteresting bits */
      if (dpd<0x16e) continue;	   /* must be canonical */
      canon=BIN2DPD[DPD2BIN[dpd]];   /* determine canonical declet */
      if (canon==dpd) continue;	   /* have canonical declet */
      /* need to replace declet */
      if (uoff>=10) {		   /* all within current word */
        encode&=~(0x3ff<<(uoff-10)); /* clear the 10 bits ready for replace */
        encode|=canon<<(uoff-10);	   /* insert the canonical form */
        DFWORD(result, inword)=encode;	/* .. and save */
        continue;
        }
      /* straddled words */
      precode=DFWORD(result, inword+1);	/* get previous */
      precode&=0xffffffff>>(10-uoff);	/* clear top bits */
      DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff)));
      encode&=0xffffffff<<uoff;		/* clear bottom bits */
      encode|=canon>>(10-uoff);		/* insert canonical */
      DFWORD(result, inword)=encode;	/* .. and save */
      } /* n */
    return result;
    } /* decCanonical */
  
  /* ------------------------------------------------------------------ */
  /* decDivide -- divide operations				      */
  /*								      */
  /*   result gets the result of dividing dfl by dfr:		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   op	    is the operation selector				      */
  /*   returns result						      */
  /*								      */
  /* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR.	      */
  /* ------------------------------------------------------------------ */
  #define DIVCOUNT  0		   /* 1 to instrument subtractions counter */
  #define DIVBASE	  BILLION	   /* the base used for divide */
  #define DIVOPLEN  DECPMAX9	   /* operand length ('digits' base 10**9) */
  #define DIVACCLEN (DIVOPLEN*3)	   /* accumulator length (ditto) */
  static decFloat * decDivide(decFloat *result, const decFloat *dfl,
  			    const decFloat *dfr, decContext *set, uInt op) {
    decFloat quotient;		   /* for remainders */
    bcdnum num;			   /* for final conversion */
    uInt	 acc[DIVACCLEN];	   /* coefficent in base-billion .. */
    uInt	 div[DIVOPLEN];		   /* divisor in base-billion .. */
    uInt	 quo[DIVOPLEN+1];	   /* quotient in base-billion .. */
    uByte	 bcdacc[(DIVOPLEN+1)*9+2]; /* for quotient in BCD, +1, +1 */
    uInt	 *msua, *msud, *msuq;	   /* -> msu of acc, div, and quo */
    Int	 divunits, accunits;	   /* lengths */
    Int	 quodigits;		   /* digits in quotient */
    uInt	 *lsua, *lsuq;		   /* -> current acc and quo lsus */
    Int	 length, multiplier;	   /* work */
    uInt	 carry, sign;		   /* .. */
    uInt	 *ua, *ud, *uq;		   /* .. */
    uByte	 *ub;			   /* .. */
    uInt	 divtop;		   /* top unit of div adjusted for estimating */
    #if DIVCOUNT
    static uInt maxcount=0;	   /* worst-seen subtractions count */
    uInt	 divcount=0;		   /* subtractions count [this divide] */
    #endif
  
    /* calculate sign */
    num.sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
  
    if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
      /* NaNs are handled as usual */
      if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
      /* one or two infinities */
      if (DFISINF(dfl)) {
        if (DFISINF(dfr)) return decInvalid(result, set); /* Two infinities bad */
        if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* as is rem */
        /* Infinity/x is infinite and quiet, even if x=0 */
        DFWORD(result, 0)=num.sign;
        return decInfinity(result, result);
        }
      /* must be x/Infinity -- remainders are lhs */
      if (op&(REMAINDER|REMNEAR)) return decCanonical(result, dfl);
      /* divides: return zero with correct sign and exponent depending */
      /* on op (Etiny for divide, 0 for divideInt) */
      decFloatZero(result);
      if (op==DIVIDEINT) DFWORD(result, 0)|=num.sign; /* add sign */
       else DFWORD(result, 0)=num.sign;	     /* zeros the exponent, too */
      return result;
      }
    /* next, handle zero operands (x/0 and 0/x) */
    if (DFISZERO(dfr)) {			     /* x/0 */
      if (DFISZERO(dfl)) {		     /* 0/0 is undefined */
        decFloatZero(result);
        DFWORD(result, 0)=DECFLOAT_qNaN;
        set->status|=DEC_Division_undefined;
        return result;
        }
      if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* bad rem */
      set->status|=DEC_Division_by_zero;
      DFWORD(result, 0)=num.sign;
      return decInfinity(result, result);	     /* x/0 -> signed Infinity */
      }
    num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr);  /* ideal exponent */
    if (DFISZERO(dfl)) {			     /* 0/x (x!=0) */
      /* if divide, result is 0 with ideal exponent; divideInt has */
      /* exponent=0, remainders give zero with lower exponent */
      if (op&DIVIDEINT) {
        decFloatZero(result);
        DFWORD(result, 0)|=num.sign;	     /* add sign */
        return result;
        }
      if (!(op&DIVIDE)) {			     /* a remainder */
        /* exponent is the minimum of the operands */
        num.exponent=MINI(GETEXPUN(dfl), GETEXPUN(dfr));
        /* if the result is zero the sign shall be sign of dfl */
        num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
        }
      bcdacc[0]=0;
      num.msd=bcdacc;			     /* -> 0 */
      num.lsd=bcdacc;			     /* .. */
      return decFinalize(result, &num, set);   /* [divide may clamp exponent] */
      } /* 0/x */
    /* [here, both operands are known to be finite and non-zero] */
  
    /* extract the operand coefficents into 'units' which are */
    /* base-billion; the lhs is high-aligned in acc and the msu of both */
    /* acc and div is at the right-hand end of array (offset length-1); */
    /* the quotient can need one more unit than the operands as digits */
    /* in it are not necessarily aligned neatly; further, the quotient */
    /* may not start accumulating until after the end of the initial */
    /* operand in acc if that is small (e.g., 1) so the accumulator */
    /* must have at least that number of units extra (at the ls end) */
    GETCOEFFBILL(dfl, acc+DIVACCLEN-DIVOPLEN);
    GETCOEFFBILL(dfr, div);
    /* zero the low uInts of acc */
    acc[0]=0;
    acc[1]=0;
    acc[2]=0;
    acc[3]=0;
    #if DOUBLE
      #if DIVOPLEN!=2
        #error Unexpected Double DIVOPLEN
      #endif
    #elif QUAD
    acc[4]=0;
    acc[5]=0;
    acc[6]=0;
    acc[7]=0;
      #if DIVOPLEN!=4
        #error Unexpected Quad DIVOPLEN
      #endif
    #endif
  
    /* set msu and lsu pointers */
    msua=acc+DIVACCLEN-1;	      /* [leading zeros removed below] */
    msuq=quo+DIVOPLEN;
    /*[loop for div will terminate because operands are non-zero] */
    for (msud=div+DIVOPLEN-1; *msud==0;) msud--;
    /* the initial least-significant unit of acc is set so acc appears */
    /* to have the same length as div. */
    /* This moves one position towards the least possible for each */
    /* iteration */
    divunits=(Int)(msud-div+1); /* precalculate */
    lsua=msua-divunits+1;	      /* initial working lsu of acc */
    lsuq=msuq;		      /* and of quo */
  
    /* set up the estimator for the multiplier; this is the msu of div, */
    /* plus two bits from the unit below (if any) rounded up by one if */
    /* there are any non-zero bits or units below that [the extra two */
    /* bits makes for a much better estimate when the top unit is small] */
    divtop=*msud<<2;
    if (divunits>1) {
      uInt *um=msud-1;
      uInt d=*um;
      if (d>=750000000) {divtop+=3; d-=750000000;}
       else if (d>=500000000) {divtop+=2; d-=500000000;}
       else if (d>=250000000) {divtop++; d-=250000000;}
      if (d) divtop++;
       else for (um--; um>=div; um--) if (*um) {
        divtop++;
        break;
        }
      } /* >1 unit */
  
    #if DECTRACE
    {Int i;
    printf("----- div=");
    for (i=divunits-1; i>=0; i--) printf("%09ld ", (LI)div[i]);
    printf("\n");}
    #endif
  
    /* now collect up to DECPMAX+1 digits in the quotient (this may */
    /* need OPLEN+1 uInts if unaligned) */
    quodigits=0;		      /* no digits yet */
    for (;; lsua--) {	      /* outer loop -- each input position */
      #if DECCHECK
      if (lsua<acc) {
        printf("Acc underrun...\n");
        break;
        }
      #endif
      #if DECTRACE
      printf("Outer: quodigits=%ld acc=", (LI)quodigits);
      for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua);
      printf("\n");
      #endif
      *lsuq=0;		      /* default unit result is 0 */
      for (;;) {		      /* inner loop -- calculate quotient unit */
        /* strip leading zero units from acc (either there initially or */
        /* from subtraction below); this may strip all if exactly 0 */
        for (; *msua==0 && msua>=lsua;) msua--;
        accunits=(Int)(msua-lsua+1);		  /* [maybe 0] */
        /* subtraction is only necessary and possible if there are as */
        /* least as many units remaining in acc for this iteration as */
        /* there are in div */
        if (accunits<divunits) {
  	if (accunits==0) msua++;		  /* restore */
  	break;
  	}
  
        /* If acc is longer than div then subtraction is definitely */
        /* possible (as msu of both is non-zero), but if they are the */
        /* same length a comparison is needed. */
        /* If a subtraction is needed then a good estimate of the */
        /* multiplier for the subtraction is also needed in order to */
        /* minimise the iterations of this inner loop because the */
        /* subtractions needed dominate division performance. */
        if (accunits==divunits) {
  	/* compare the high divunits of acc and div: */
  	/* acc<div:  this quotient unit is unchanged; subtraction */
  	/*	     will be possible on the next iteration */
  	/* acc==div: quotient gains 1, set acc=0 */
  	/* acc>div:  subtraction necessary at this position */
  	for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break;
  	/* [now at first mismatch or lsu] */
  	if (*ud>*ua) break;			  /* next time... */
  	if (*ud==*ua) {				  /* all compared equal */
  	  *lsuq+=1;				  /* increment result */
  	  msua=lsua;				  /* collapse acc units */
  	  *msua=0;				  /* .. to a zero */
  	  break;
  	  }
  
  	/* subtraction necessary; estimate multiplier [see above] */
  	/* if both *msud and *msua are small it is cost-effective to */
  	/* bring in part of the following units (if any) to get a */
  	/* better estimate (assume some other non-zero in div) */
  	#define DIVLO 1000000U
  	#define DIVHI (DIVBASE/DIVLO)
  	#if DECUSE64
  	  if (divunits>1) {
  	    /* there cannot be a *(msud-2) for DECDOUBLE so next is */
  	    /* an exact calculation unless DECQUAD (which needs to */
  	    /* assume bits out there if divunits>2) */
  	    uLong mul=(uLong)*msua * DIVBASE + *(msua-1);
  	    uLong div=(uLong)*msud * DIVBASE + *(msud-1);
  	    #if QUAD
  	    if (divunits>2) div++;
  	    #endif
  	    mul/=div;
  	    multiplier=(Int)mul;
  	    }
  	   else multiplier=*msua/(*msud);
  	#else
  	  if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
  	    multiplier=(*msua*DIVHI + *(msua-1)/DIVLO)
  		      /(*msud*DIVHI + *(msud-1)/DIVLO +1);
  	    }
  	   else multiplier=(*msua<<2)/divtop;
  	#endif
  	}
         else {					  /* accunits>divunits */
  	/* msud is one unit 'lower' than msua, so estimate differently */
  	#if DECUSE64
  	  uLong mul;
  	  /* as before, bring in extra digits if possible */
  	  if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
  	    mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI
  	       + *(msua-2)/DIVLO;
  	    mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1);
  	    }
  	   else if (divunits==1) {
  	    mul=(uLong)*msua * DIVBASE + *(msua-1);
  	    mul/=*msud;	      /* no more to the right */
  	    }
  	   else {
  	    mul=(uLong)(*msua) * (uInt)(DIVBASE<<2) + (*(msua-1)<<2);
  	    mul/=divtop;      /* [divtop already allows for sticky bits] */
  	    }
  	  multiplier=(Int)mul;
  	#else
  	  multiplier=*msua * ((DIVBASE<<2)/divtop);
  	#endif
  	}
        if (multiplier==0) multiplier=1;		  /* marginal case */
        *lsuq+=multiplier;
  
        #if DIVCOUNT
        /* printf("Multiplier: %ld\n", (LI)multiplier); */
        divcount++;
        #endif
  
        /* Carry out the subtraction  acc-(div*multiplier); for each */
        /* unit in div, do the multiply, split to units (see */
        /* decFloatMultiply for the algorithm), and subtract from acc */
        #define DIVMAGIC	2305843009U		  /* 2**61/10**9 */
        #define DIVSHIFTA 29
        #define DIVSHIFTB 32
        carry=0;
        for (ud=div, ua=lsua; ud<=msud; ud++, ua++) {
  	uInt lo, hop;
  	#if DECUSE64
  	  uLong sub=(uLong)multiplier*(*ud)+carry;
  	  if (sub<DIVBASE) {
  	    carry=0;
  	    lo=(uInt)sub;
  	    }
  	   else {
  	    hop=(uInt)(sub>>DIVSHIFTA);
  	    carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB);
  	    /* the estimate is now in hi; now calculate sub-hi*10**9 */
  	    /* to get the remainder (which will be <DIVBASE)) */
  	    lo=(uInt)sub;
  	    lo-=carry*DIVBASE;			  /* low word of result */
  	    if (lo>=DIVBASE) {
  	      lo-=DIVBASE;			  /* correct by +1 */
  	      carry++;
  	      }
  	    }
  	#else /* 32-bit */
  	  uInt hi;
  	  /* calculate multiplier*(*ud) into hi and lo */
  	  LONGMUL32HI(hi, *ud, multiplier);	  /* get the high word */
  	  lo=multiplier*(*ud);			  /* .. and the low */
  	  lo+=carry;				  /* add the old hi */
  	  carry=hi+(lo<carry);			  /* .. with any carry */
  	  if (carry || lo>=DIVBASE) {		  /* split is needed */
  	    hop=(carry<<3)+(lo>>DIVSHIFTA);	  /* hi:lo/2**29 */
  	    LONGMUL32HI(carry, hop, DIVMAGIC);	  /* only need the high word */
  	    /* [DIVSHIFTB is 32, so carry can be used directly] */
  	    /* the estimate is now in carry; now calculate hi:lo-est*10**9; */
  	    /* happily the top word of the result is irrelevant because it */
  	    /* will always be zero so this needs only one multiplication */
  	    lo-=(carry*DIVBASE);
  	    /* the correction here will be at most +1; do it */
  	    if (lo>=DIVBASE) {
  	      lo-=DIVBASE;
  	      carry++;
  	      }
  	    }
  	#endif
  	if (lo>*ua) {		   /* borrow needed */
  	  *ua+=DIVBASE;
  	  carry++;
  	  }
  	*ua-=lo;
  	} /* ud loop */
        if (carry) *ua-=carry;	   /* accdigits>divdigits [cannot borrow] */
        } /* inner loop */
  
      /* the outer loop terminates when there is either an exact result */
      /* or enough digits; first update the quotient digit count and */
      /* pointer (if any significant digits) */
      #if DECTRACE
      if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq);
      #endif
      if (quodigits) {
        quodigits+=9;		   /* had leading unit earlier */
        lsuq--;
        if (quodigits>DECPMAX+1) break;	/* have enough */
        }
       else if (*lsuq) {		   /* first quotient digits */
        const uInt *pow;
        for (pow=DECPOWERS; *lsuq>=*pow; pow++) quodigits++;
        lsuq--;
        /* [cannot have >DECPMAX+1 on first unit] */
        }
  
      if (*msua!=0) continue;	   /* not an exact result */
      /* acc is zero iff used all of original units and zero down to lsua */
      /* (must also continue to original lsu for correct quotient length) */
      if (lsua>acc+DIVACCLEN-DIVOPLEN) continue;
      for (; msua>lsua && *msua==0;) msua--;
      if (*msua==0 && msua==lsua) break;
      } /* outer loop */
  
    /* all of the original operand in acc has been covered at this point */
    /* quotient now has at least DECPMAX+2 digits */
    /* *msua is now non-0 if inexact and sticky bits */
    /* lsuq is one below the last uint of the quotient */
    lsuq++;			   /* set -> true lsu of quo */
    if (*msua) *lsuq|=1;		   /* apply sticky bit */
  
    /* quo now holds the (unrounded) quotient in base-billion; one */
    /* base-billion 'digit' per uInt. */
    #if DECTRACE
    printf("DivQuo:");
    for (uq=msuq; uq>=lsuq; uq--) printf(" %09ld", (LI)*uq);
    printf("\n");
    #endif
  
    /* Now convert to BCD for rounding and cleanup, starting from the */
    /* most significant end [offset by one into bcdacc to leave room */
    /* for a possible carry digit if rounding for REMNEAR is needed] */
    for (uq=msuq, ub=bcdacc+1; uq>=lsuq; uq--, ub+=9) {
      uInt top, mid, rem;			/* work */
      if (*uq==0) {			/* no split needed */
        UINTAT(ub)=0;			/* clear 9 BCD8s */
        UINTAT(ub+4)=0;			/* .. */
        *(ub+8)=0;			/* .. */
        continue;
        }
      /* *uq is non-zero -- split the base-billion digit into */
      /* hi, mid, and low three-digits */
      #define divsplit9 1000000		/* divisor */
      #define divsplit6 1000		/* divisor */
      /* The splitting is done by simple divides and remainders, */
      /* assuming the compiler will optimize these [GCC does] */
      top=*uq/divsplit9;
      rem=*uq%divsplit9;
      mid=rem/divsplit6;
      rem=rem%divsplit6;
      /* lay out the nine BCD digits (plus one unwanted byte) */
      UINTAT(ub)	=UINTAT(&BIN2BCD8[top*4]);
      UINTAT(ub+3)=UINTAT(&BIN2BCD8[mid*4]);
      UINTAT(ub+6)=UINTAT(&BIN2BCD8[rem*4]);
      } /* BCD conversion loop */
    ub--;					/* -> lsu */
  
    /* complete the bcdnum; quodigits is correct, so the position of */
    /* the first non-zero is known */
    num.msd=bcdacc+1+(msuq-lsuq+1)*9-quodigits;
    num.lsd=ub;
  
    /* make exponent adjustments, etc */
    if (lsua<acc+DIVACCLEN-DIVOPLEN) {	/* used extra digits */
      num.exponent-=(Int)((acc+DIVACCLEN-DIVOPLEN-lsua)*9);
      /* if the result was exact then there may be up to 8 extra */
      /* trailing zeros in the overflowed quotient final unit */
      if (*msua==0) {
        for (; *ub==0;) ub--;		/* drop zeros */
        num.exponent+=(Int)(num.lsd-ub);	/* and adjust exponent */
        num.lsd=ub;
        }
      } /* adjustment needed */
  
    #if DIVCOUNT
    if (divcount>maxcount) {		/* new high-water nark */
      maxcount=divcount;
      printf("DivNewMaxCount: %ld\n", (LI)maxcount);
      }
    #endif
  
    if (op&DIVIDE) return decFinalize(result, &num, set); /* all done */
  
    /* Is DIVIDEINT or a remainder; there is more to do -- first form */
    /* the integer (this is done 'after the fact', unlike as in */
    /* decNumber, so as not to tax DIVIDE) */
  
    /* The first non-zero digit will be in the first 9 digits, known */
    /* from quodigits and num.msd, so there is always space for DECPMAX */
    /* digits */
  
    length=(Int)(num.lsd-num.msd+1);
    /*printf("Length exp: %ld %ld\n", (LI)length, (LI)num.exponent); */
  
    if (length+num.exponent>DECPMAX) { /* cannot fit */
      decFloatZero(result);
      DFWORD(result, 0)=DECFLOAT_qNaN;
      set->status|=DEC_Division_impossible;
      return result;
      }
  
    if (num.exponent>=0) {	   /* already an int, or need pad zeros */
      for (ub=num.lsd+1; ub<=num.lsd+num.exponent; ub++) *ub=0;
      num.lsd+=num.exponent;
      }
     else {			   /* too long: round or truncate needed */
      Int drop=-num.exponent;
      if (!(op&REMNEAR)) {	   /* simple truncate */
        num.lsd-=drop;
        if (num.lsd<num.msd) {	   /* truncated all */
  	num.lsd=num.msd;	   /* make 0 */
  	*num.lsd=0;		   /* .. [sign still relevant] */
  	}
        }
       else {			   /* round to nearest even [sigh] */
        /* round-to-nearest, in-place; msd is at or to right of bcdacc+1 */
        /* (this is a special case of Quantize -- q.v. for commentary) */
        uByte *roundat;		   /* -> re-round digit */
        uByte reround;		   /* reround value */
        *(num.msd-1)=0;		   /* in case of left carry, or make 0 */
        if (drop<length) roundat=num.lsd-drop+1;
         else if (drop==length) roundat=num.msd;
         else roundat=num.msd-1;	   /* [-> 0] */
        reround=*roundat;
        for (ub=roundat+1; ub<=num.lsd; ub++) {
  	if (*ub!=0) {
  	  reround=DECSTICKYTAB[reround];
  	  break;
  	  }
  	} /* check stickies */
        if (roundat>num.msd) num.lsd=roundat-1;
         else {
  	num.msd--;			     /* use the 0 .. */
  	num.lsd=num.msd;		     /* .. at the new MSD place */
  	}
        if (reround!=0) {			     /* discarding non-zero */
  	uInt bump=0;
  	/* rounding is DEC_ROUND_HALF_EVEN always */
  	if (reround>5) bump=1;		     /* >0.5 goes up */
  	 else if (reround==5)		     /* exactly 0.5000 .. */
  	  bump=*(num.lsd) & 0x01;	     /* .. up iff [new] lsd is odd */
  	if (bump!=0) {			     /* need increment */
  	  /* increment the coefficient; this might end up with 1000... */
  	  ub=num.lsd;
  	  for (; UINTAT(ub-3)==0x09090909; ub-=4) UINTAT(ub-3)=0;
  	  for (; *ub==9; ub--) *ub=0;	     /* at most 3 more */
  	  *ub+=1;
  	  if (ub<num.msd) num.msd--;	     /* carried */
  	  } /* bump needed */
  	} /* reround!=0 */
        } /* remnear */
      } /* round or truncate needed */
    num.exponent=0;			     /* all paths */
    /*decShowNum(&num, "int"); */
  
    if (op&DIVIDEINT) return decFinalize(result, &num, set); /* all done */
  
    /* Have a remainder to calculate */
    decFinalize(&quotient, &num, set);	     /* lay out the integer so far */
    DFWORD(&quotient, 0)^=DECFLOAT_Sign;	     /* negate it */
    sign=DFWORD(dfl, 0);			     /* save sign of dfl */
    decFloatFMA(result, &quotient, dfr, dfl, set);
    if (!DFISZERO(result)) return result;
    /* if the result is zero the sign shall be sign of dfl */
    DFWORD(&quotient, 0)=sign;		     /* construct decFloat of sign */
    return decFloatCopySign(result, result, &quotient);
    } /* decDivide */
  
  /* ------------------------------------------------------------------ */
  /* decFiniteMultiply -- multiply two finite decFloats		      */
  /*								      */
  /*   num    gets the result of multiplying dfl and dfr		      */
  /*   bcdacc .. with the coefficient in this array		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*								      */
  /* This effects the multiplication of two decFloats, both known to be */
  /* finite, leaving the result in a bcdnum ready for decFinalize (for  */
  /* use in Multiply) or in a following addition (FMA).		      */
  /*								      */
  /* bcdacc must have space for at least DECPMAX9*18+1 bytes.	      */
  /* No error is possible and no status is set.			      */
  /* ------------------------------------------------------------------ */
  /* This routine has two separate implementations of the core */
  /* multiplication; both using base-billion.  One uses only 32-bit */
  /* variables (Ints and uInts) or smaller; the other uses uLongs (for */
  /* multiplication and addition only).  Both implementations cover */
  /* both arithmetic sizes (DOUBLE and QUAD) in order to allow timing */
  /* comparisons.	 In any one compilation only one implementation for */
  /* each size can be used, and if DECUSE64 is 0 then use of the 32-bit */
  /* version is forced. */
  /* */
  /* Historical note: an earlier version of this code also supported the */
  /* 256-bit format and has been preserved.  That is somewhat trickier */
  /* during lazy carry splitting because the initial quotient estimate */
  /* (est) can exceed 32 bits. */
  
  #define MULTBASE  BILLION	   /* the base used for multiply */
  #define MULOPLEN  DECPMAX9	   /* operand length ('digits' base 10**9) */
  #define MULACCLEN (MULOPLEN*2)		    /* accumulator length (ditto) */
  #define LEADZEROS (MULACCLEN*9 - DECPMAX*2) /* leading zeros always */
  
  /* Assertions: exponent not too large and MULACCLEN is a multiple of 4 */
  #if DECEMAXD>9
    #error Exponent may overflow when doubled for Multiply
  #endif
  #if MULACCLEN!=(MULACCLEN/4)*4
    /* This assumption is used below only for initialization */
    #error MULACCLEN is not a multiple of 4
  #endif
  
  static void decFiniteMultiply(bcdnum *num, uByte *bcdacc,
  			      const decFloat *dfl, const decFloat *dfr) {
    uInt	 bufl[MULOPLEN];	   /* left  coefficient (base-billion) */
    uInt	 bufr[MULOPLEN];	   /* right coefficient (base-billion) */
    uInt	 *ui, *uj;		   /* work */
    uByte	 *ub;			   /* .. */
  
    #if DECUSE64
    uLong	 accl[MULACCLEN];	   /* lazy accumulator (base-billion+) */
    uLong	 *pl;			   /* work -> lazy accumulator */
    uInt	 acc[MULACCLEN];	   /* coefficent in base-billion .. */
    #else
    uInt	 acc[MULACCLEN*2];	   /* accumulator in base-billion .. */
    #endif
    uInt	 *pa;			   /* work -> accumulator */
    /*printf("Base10**9: OpLen=%d MulAcclen=%d\n", OPLEN, MULACCLEN); */
  
    /* Calculate sign and exponent */
    num->sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
    num->exponent=GETEXPUN(dfl)+GETEXPUN(dfr); /* [see assertion above] */
  
    /* Extract the coefficients and prepare the accumulator */
    /* the coefficients of the operands are decoded into base-billion */
    /* numbers in uInt arrays (bufl and bufr, LSD at offset 0) of the */
    /* appropriate size. */
    GETCOEFFBILL(dfl, bufl);
    GETCOEFFBILL(dfr, bufr);
    #if DECTRACE && 0
      printf("CoeffbL:");
      for (ui=bufl+MULOPLEN-1; ui>=bufl; ui--) printf(" %08lx", (LI)*ui);
      printf("\n");
      printf("CoeffbR:");
      for (uj=bufr+MULOPLEN-1; uj>=bufr; uj--) printf(" %08lx", (LI)*uj);
      printf("\n");
    #endif
  
    /* start the 64-bit/32-bit differing paths... */
  #if DECUSE64
  
    /* zero the accumulator */
    #if MULACCLEN==4
      accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0;
    #else					     /* use a loop */
      /* MULACCLEN is a multiple of four, asserted above */
      for (pl=accl; pl<accl+MULACCLEN; pl+=4) {
        *pl=0; *(pl+1)=0; *(pl+2)=0; *(pl+3)=0;/* [reduce overhead] */
        } /* pl */
    #endif
  
    /* Effect the multiplication */
    /* The multiplcation proceeds using MFC's lazy-carry resolution */
    /* algorithm from decNumber.	First, the multiplication is */
    /* effected, allowing accumulation of the partial products (which */
    /* are in base-billion at each column position) into 64 bits */
    /* without resolving back to base=billion after each addition. */
    /* These 64-bit numbers (which may contain up to 19 decimal digits) */
    /* are then split using the Clark & Cowlishaw algorithm (see below). */
    /* [Testing for 0 in the inner loop is not really a 'win'] */
    for (ui=bufr; ui<bufr+MULOPLEN; ui++) { /* over each item in rhs */
      if (*ui==0) continue;		  /* product cannot affect result */
      pl=accl+(ui-bufr);			  /* where to add the lhs */
      for (uj=bufl; uj<bufl+MULOPLEN; uj++, pl++) { /* over each item in lhs */
        /* if (*uj==0) continue;		  // product cannot affect result */
        *pl+=((uLong)*ui)*(*uj);
        } /* uj */
      } /* ui */
  
    /* The 64-bit carries must now be resolved; this means that a */
    /* quotient/remainder has to be calculated for base-billion (1E+9). */
    /* For this, Clark & Cowlishaw's quotient estimation approach (also */
    /* used in decNumber) is needed, because 64-bit divide is generally */
    /* extremely slow on 32-bit machines, and may be slower than this */
    /* approach even on 64-bit machines.	This algorithm splits X */
    /* using: */
    /* */
    /*   magic=2**(A+B)/1E+9;   // 'magic number' */
    /*   hop=X/2**A;	      // high order part of X (by shift) */
    /*   est=magic*hop/2**B     // quotient estimate (may be low by 1) */
    /* */
    /* A and B are quite constrained; hop and magic must fit in 32 bits, */
    /* and 2**(A+B) must be as large as possible (which is 2**61 if */
    /* magic is to fit).	Further, maxX increases with the length of */
    /* the operands (and hence the number of partial products */
    /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
    /* */
    /* It can be shown that when OPLEN is 2 then the maximum error in */
    /* the estimated quotient is <1, but for larger maximum x the */
    /* maximum error is above 1 so a correction that is >1 may be */
    /* needed.  Values of A and B are chosen to satisfy the constraints */
    /* just mentioned while minimizing the maximum error (and hence the */
    /* maximum correction), as shown in the following table: */
    /* */
    /*   Type    OPLEN   A   B	 maxX	 maxError  maxCorrection */
    /*   --------------------------------------------------------- */
    /*   DOUBLE	 2    29  32  <2*10**18	   0.63	      1 */
    /*   QUAD	 4    30  31  <4*10**18	   1.17	      2 */
    /* */
    /* In the OPLEN==2 case there is most choice, but the value for B */
    /* of 32 has a big advantage as then the calculation of the */
    /* estimate requires no shifting; the compiler can extract the high */
    /* word directly after multiplying magic*hop. */
    #define MULMAGIC 2305843009U		/* 2**61/10**9	[both cases] */
    #if DOUBLE
      #define MULSHIFTA 29
      #define MULSHIFTB 32
    #elif QUAD
      #define MULSHIFTA 30
      #define MULSHIFTB 31
    #else
      #error Unexpected type
    #endif
  
    #if DECTRACE
    printf("MulAccl:");
    for (pl=accl+MULACCLEN-1; pl>=accl; pl--)
      printf(" %08lx:%08lx", (LI)(*pl>>32), (LI)(*pl&0xffffffff));
    printf("\n");
    #endif
  
    for (pl=accl, pa=acc; pl<accl+MULACCLEN; pl++, pa++) { /* each column position */
      uInt lo, hop;			/* work */
      uInt est;				/* cannot exceed 4E+9 */
      if (*pl>MULTBASE) {
        /* *pl holds a binary number which needs to be split */
        hop=(uInt)(*pl>>MULSHIFTA);
        est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB);
        /* the estimate is now in est; now calculate hi:lo-est*10**9; */
        /* happily the top word of the result is irrelevant because it */
        /* will always be zero so this needs only one multiplication */
        lo=(uInt)(*pl-((uLong)est*MULTBASE));  /* low word of result */
        /* If QUAD, the correction here could be +2 */
        if (lo>=MULTBASE) {
  	lo-=MULTBASE;			/* correct by +1 */
  	est++;
  	#if QUAD
  	/* may need to correct by +2 */
  	if (lo>=MULTBASE) {
  	  lo-=MULTBASE;
  	  est++;
  	  }
  	#endif
  	}
        /* finally place lo as the new coefficient 'digit' and add est to */
        /* the next place up [this is safe because this path is never */
        /* taken on the final iteration as *pl will fit] */
        *pa=lo;
        *(pl+1)+=est;
        } /* *pl needed split */
       else {				/* *pl<MULTBASE */
        *pa=(uInt)*pl;			/* just copy across */
        }
      } /* pl loop */
  
  #else  /* 32-bit */
    for (pa=acc;; pa+=4) {		     /* zero the accumulator */
      *pa=0; *(pa+1)=0; *(pa+2)=0; *(pa+3)=0;  /* [reduce overhead] */
      if (pa==acc+MULACCLEN*2-4) break;	     /* multiple of 4 asserted */
      } /* pa */
  
    /* Effect the multiplication */
    /* uLongs are not available (and in particular, there is no uLong */
    /* divide) but it is still possible to use MFC's lazy-carry */
    /* resolution algorithm from decNumber.  First, the multiplication */
    /* is effected, allowing accumulation of the partial products */
    /* (which are in base-billion at each column position) into 64 bits */
    /* [with the high-order 32 bits in each position being held at */
    /* offset +ACCLEN from the low-order 32 bits in the accumulator]. */
    /* These 64-bit numbers (which may contain up to 19 decimal digits) */
    /* are then split using the Clark & Cowlishaw algorithm (see */
    /* below). */
    for (ui=bufr;; ui++) {		/* over each item in rhs */
      uInt hi, lo;			/* words of exact multiply result */
      pa=acc+(ui-bufr);			/* where to add the lhs */
      for (uj=bufl;; uj++, pa++) {	/* over each item in lhs */
        LONGMUL32HI(hi, *ui, *uj);	/* calculate product of digits */
        lo=(*ui)*(*uj);			/* .. */
        *pa+=lo;				/* accumulate low bits and .. */
        *(pa+MULACCLEN)+=hi+(*pa<lo);	/* .. high bits with any carry */
        if (uj==bufl+MULOPLEN-1) break;
        }
      if (ui==bufr+MULOPLEN-1) break;
      }
  
    /* The 64-bit carries must now be resolved; this means that a */
    /* quotient/remainder has to be calculated for base-billion (1E+9). */
    /* For this, Clark & Cowlishaw's quotient estimation approach (also */
    /* used in decNumber) is needed, because 64-bit divide is generally */
    /* extremely slow on 32-bit machines.	 This algorithm splits X */
    /* using: */
    /* */
    /*   magic=2**(A+B)/1E+9;   // 'magic number' */
    /*   hop=X/2**A;	      // high order part of X (by shift) */
    /*   est=magic*hop/2**B     // quotient estimate (may be low by 1) */
    /* */
    /* A and B are quite constrained; hop and magic must fit in 32 bits, */
    /* and 2**(A+B) must be as large as possible (which is 2**61 if */
    /* magic is to fit).	Further, maxX increases with the length of */
    /* the operands (and hence the number of partial products */
    /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
    /* */
    /* It can be shown that when OPLEN is 2 then the maximum error in */
    /* the estimated quotient is <1, but for larger maximum x the */
    /* maximum error is above 1 so a correction that is >1 may be */
    /* needed.  Values of A and B are chosen to satisfy the constraints */
    /* just mentioned while minimizing the maximum error (and hence the */
    /* maximum correction), as shown in the following table: */
    /* */
    /*   Type    OPLEN   A   B	 maxX	 maxError  maxCorrection */
    /*   --------------------------------------------------------- */
    /*   DOUBLE	 2    29  32  <2*10**18	   0.63	      1 */
    /*   QUAD	 4    30  31  <4*10**18	   1.17	      2 */
    /* */
    /* In the OPLEN==2 case there is most choice, but the value for B */
    /* of 32 has a big advantage as then the calculation of the */
    /* estimate requires no shifting; the high word is simply */
    /* calculated from multiplying magic*hop. */
    #define MULMAGIC 2305843009U		/* 2**61/10**9	[both cases] */
    #if DOUBLE
      #define MULSHIFTA 29
      #define MULSHIFTB 32
    #elif QUAD
      #define MULSHIFTA 30
      #define MULSHIFTB 31
    #else
      #error Unexpected type
    #endif
  
    #if DECTRACE
    printf("MulHiLo:");
    for (pa=acc+MULACCLEN-1; pa>=acc; pa--)
      printf(" %08lx:%08lx", (LI)*(pa+MULACCLEN), (LI)*pa);
    printf("\n");
    #endif
  
    for (pa=acc;; pa++) {			/* each low uInt */
      uInt hi, lo;			/* words of exact multiply result */
      uInt hop, estlo;			/* work */
      #if QUAD
      uInt esthi;				/* .. */
      #endif
  
      lo=*pa;
      hi=*(pa+MULACCLEN);			/* top 32 bits */
      /* hi and lo now hold a binary number which needs to be split */
  
      #if DOUBLE
        hop=(hi<<3)+(lo>>MULSHIFTA);	/* hi:lo/2**29 */
        LONGMUL32HI(estlo, hop, MULMAGIC);/* only need the high word */
        /* [MULSHIFTB is 32, so estlo can be used directly] */
        /* the estimate is now in estlo; now calculate hi:lo-est*10**9; */
        /* happily the top word of the result is irrelevant because it */
        /* will always be zero so this needs only one multiplication */
        lo-=(estlo*MULTBASE);
        /* esthi=0;			// high word is ignored below */
        /* the correction here will be at most +1; do it */
        if (lo>=MULTBASE) {
  	lo-=MULTBASE;
  	estlo++;
  	}
      #elif QUAD
        hop=(hi<<2)+(lo>>MULSHIFTA);	/* hi:lo/2**30 */
        LONGMUL32HI(esthi, hop, MULMAGIC);/* shift will be 31 .. */
        estlo=hop*MULMAGIC;		/* .. so low word needed */
        estlo=(esthi<<1)+(estlo>>MULSHIFTB); /* [just the top bit] */
        /* esthi=0;			// high word is ignored below */
        lo-=(estlo*MULTBASE);		/* as above */
        /* the correction here could be +1 or +2 */
        if (lo>=MULTBASE) {
  	lo-=MULTBASE;
  	estlo++;
  	}
        if (lo>=MULTBASE) {
  	lo-=MULTBASE;
  	estlo++;
  	}
      #else
        #error Unexpected type
      #endif
  
      /* finally place lo as the new accumulator digit and add est to */
      /* the next place up; this latter add could cause a carry of 1 */
      /* to the high word of the next place */
      *pa=lo;
      *(pa+1)+=estlo;
      /* esthi is always 0 for DOUBLE and QUAD so this is skipped */
      /* *(pa+1+MULACCLEN)+=esthi; */
      if (*(pa+1)<estlo) *(pa+1+MULACCLEN)+=1; /* carry */
      if (pa==acc+MULACCLEN-2) break;	     /* [MULACCLEN-1 will never need split] */
      } /* pa loop */
  #endif
  
    /* At this point, whether using the 64-bit or the 32-bit paths, the */
    /* accumulator now holds the (unrounded) result in base-billion; */
    /* one base-billion 'digit' per uInt. */
    #if DECTRACE
    printf("MultAcc:");
    for (pa=acc+MULACCLEN-1; pa>=acc; pa--) printf(" %09ld", (LI)*pa);
    printf("\n");
    #endif
  
    /* Now convert to BCD for rounding and cleanup, starting from the */
    /* most significant end */
    pa=acc+MULACCLEN-1;
    if (*pa!=0) num->msd=bcdacc+LEADZEROS;/* drop known lead zeros */
     else {				/* >=1 word of leading zeros */
      num->msd=bcdacc;			/* known leading zeros are gone */
      pa--;				/* skip first word .. */
      for (; *pa==0; pa--) if (pa==acc) break; /* .. and any more leading 0s */
      }
    for (ub=bcdacc;; pa--, ub+=9) {
      if (*pa!=0) {			/* split(s) needed */
        uInt top, mid, rem;		/* work */
        /* *pa is non-zero -- split the base-billion acc digit into */
        /* hi, mid, and low three-digits */
        #define mulsplit9 1000000		/* divisor */
        #define mulsplit6 1000		/* divisor */
        /* The splitting is done by simple divides and remainders, */
        /* assuming the compiler will optimize these where useful */
        /* [GCC does] */
        top=*pa/mulsplit9;
        rem=*pa%mulsplit9;
        mid=rem/mulsplit6;
        rem=rem%mulsplit6;
        /* lay out the nine BCD digits (plus one unwanted byte) */
        UINTAT(ub)  =UINTAT(&BIN2BCD8[top*4]);
        UINTAT(ub+3)=UINTAT(&BIN2BCD8[mid*4]);
        UINTAT(ub+6)=UINTAT(&BIN2BCD8[rem*4]);
        }
       else {				/* *pa==0 */
        UINTAT(ub)=0;			/* clear 9 BCD8s */
        UINTAT(ub+4)=0;			/* .. */
        *(ub+8)=0;			/* .. */
        }
      if (pa==acc) break;
      } /* BCD conversion loop */
  
    num->lsd=ub+8;			/* complete the bcdnum .. */
  
    #if DECTRACE
    decShowNum(num, "postmult");
    decFloatShow(dfl, "dfl");
    decFloatShow(dfr, "dfr");
    #endif
    return;
    } /* decFiniteMultiply */
  
  /* ------------------------------------------------------------------ */
  /* decFloatAbs -- absolute value, heeding NaNs, etc.		      */
  /*								      */
  /*   result gets the canonicalized df with sign 0		      */
  /*   df	    is the decFloat to abs				      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* This has the same effect as decFloatPlus unless df is negative,    */
  /* in which case it has the same effect as decFloatMinus.  The	      */
  /* effect is also the same as decFloatCopyAbs except that NaNs are    */
  /* handled normally (the sign of a NaN is not affected, and an sNaN   */
  /* will signal) and the result will be canonical.		      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatAbs(decFloat *result, const decFloat *df,
  		       decContext *set) {
    if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
    decCanonical(result, df);		/* copy and check */
    DFBYTE(result, 0)&=~0x80;		/* zero sign bit */
    return result;
    } /* decFloatAbs */
  
  /* ------------------------------------------------------------------ */
  /* decFloatAdd -- add two decFloats				      */
  /*								      */
  /*   result gets the result of adding dfl and dfr:		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatAdd(decFloat *result,
  		       const decFloat *dfl, const decFloat *dfr,
  		       decContext *set) {
    bcdnum num;			   /* for final conversion */
    Int	 expl, expr;		   /* left and right exponents */
    uInt	 *ui, *uj;		   /* work */
    uByte	 *ub;			   /* .. */
  
    uInt sourhil, sourhir;	   /* top words from source decFloats */
  				   /* [valid only until specials */
  				   /* handled or exponents decoded] */
    uInt diffsign;		   /* non-zero if signs differ */
    uInt carry;			   /* carry: 0 or 1 before add loop */
    Int  overlap;			   /* coefficient overlap (if full) */
    /* the following buffers hold coefficients with various alignments */
    /* (see commentary and diagrams below) */
    uByte acc[4+2+DECPMAX*3+8];
    uByte buf[4+2+DECPMAX*2];
    uByte *umsd, *ulsd;		   /* local MSD and LSD pointers */
  
    #if DECLITEND
      #define CARRYPAT 0x01000000	   /* carry=1 pattern */
    #else
      #define CARRYPAT 0x00000001	   /* carry=1 pattern */
    #endif
  
    /* Start decoding the arguments */
    /* the initial exponents are placed into the opposite Ints to */
    /* that which might be expected; there are two sets of data to */
    /* keep track of (each decFloat and the corresponding exponent), */
    /* and this scheme means that at the swap point (after comparing */
    /* exponents) only one pair of words needs to be swapped */
    /* whichever path is taken (thereby minimising worst-case path) */
    sourhil=DFWORD(dfl, 0);	   /* LHS top word */
    expr=DECCOMBEXP[sourhil>>26];	   /* get exponent high bits (in place) */
    sourhir=DFWORD(dfr, 0);	   /* RHS top word */
    expl=DECCOMBEXP[sourhir>>26];
  
    diffsign=(sourhil^sourhir)&DECFLOAT_Sign;
  
    if (EXPISSPECIAL(expl | expr)) { /* either is special? */
      if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
      /* one or two infinities */
      /* two infinities with different signs is invalid */
      if (diffsign && DFISINF(dfl) && DFISINF(dfr))
        return decInvalid(result, set);
      if (DFISINF(dfl)) return decInfinity(result, dfl); /* LHS is infinite */
      return decInfinity(result, dfr);		       /* RHS must be Infinite */
      }
  
    /* Here when both arguments are finite */
  
    /* complete exponent gathering (keeping swapped) */
    expr+=GETECON(dfl)-DECBIAS;	   /* .. + continuation and unbias */
    expl+=GETECON(dfr)-DECBIAS;
    /* here expr has exponent from lhs, and vice versa */
  
    /* now swap either exponents or argument pointers */
    if (expl<=expr) {
      /* original left is bigger */
      Int expswap=expl;
      expl=expr;
      expr=expswap;
      /* printf("left bigger\n"); */
      }
     else {
      const decFloat *dfswap=dfl;
      dfl=dfr;
      dfr=dfswap;
      /* printf("right bigger\n"); */
      }
    /* [here dfl and expl refer to the datum with the larger exponent, */
    /* of if the exponents are equal then the original LHS argument] */
  
    /* if lhs is zero then result will be the rhs (now known to have */
    /* the smaller exponent), which also may need to be tested for zero */
    /* for the weird IEEE 754 sign rules */
    if (DFISZERO(dfl)) {
      decCanonical(result, dfr);		     /* clean copy */
      /* "When the sum of two operands with opposite signs is */
      /* exactly zero, the sign of that sum shall be '+' in all */
      /* rounding modes except round toward -Infinity, in which */
      /* mode that sign shall be '-'." */
      if (diffsign && DFISZERO(result)) {
        DFWORD(result, 0)&=~DECFLOAT_Sign;     /* assume sign 0 */
        if (set->round==DEC_ROUND_FLOOR) DFWORD(result, 0)|=DECFLOAT_Sign;
        }
      return result;
      } /* numfl is zero */
    /* [here, LHS is non-zero; code below assumes that] */
  
    /* Coefficients layout during the calculations to follow: */
    /* */
    /*	   Overlap case: */
    /*	   +------------------------------------------------+ */
    /* acc:  |0000|      coeffa	   | tail B |		    | */
    /*	   +------------------------------------------------+ */
    /* buf:  |0000| pad0s |      coeffb	    |		    | */
    /*	   +------------------------------------------------+ */
    /* */
    /*	   Touching coefficients or gap: */
    /*	   +------------------------------------------------+ */
    /* acc:  |0000|      coeffa	   | gap |	coeffb	    | */
    /*	   +------------------------------------------------+ */
    /*	   [buf not used or needed; gap clamped to Pmax] */
  
    /* lay out lhs coefficient into accumulator; this starts at acc+4 */
    /* for decDouble or acc+6 for decQuad so the LSD is word- */
    /* aligned; the top word gap is there only in case a carry digit */
    /* is prefixed after the add -- it does not need to be zeroed */
    #if DOUBLE
      #define COFF 4			/* offset into acc */
    #elif QUAD
      USHORTAT(acc+4)=0;			/* prefix 00 */
      #define COFF 6			/* offset into acc */
    #endif
  
    GETCOEFF(dfl, acc+COFF);		/* decode from decFloat */
    ulsd=acc+COFF+DECPMAX-1;
    umsd=acc+4;				/* [having this here avoids */
  					/* weird GCC optimizer failure] */
    #if DECTRACE
    {bcdnum tum;
    tum.msd=umsd;
    tum.lsd=ulsd;
    tum.exponent=expl;
    tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
    decShowNum(&tum, "dflx");}
    #endif
  
    /* if signs differ, take ten's complement of lhs (here the */
    /* coefficient is subtracted from all-nines; the 1 is added during */
    /* the later add cycle -- zeros to the right do not matter because */
    /* the complement of zero is zero); these are fixed-length inverts */
    /* where the lsd is known to be at a 4-byte boundary (so no borrow */
    /* possible) */
    carry=0;				/* assume no carry */
    if (diffsign) {
      carry=CARRYPAT;			/* for +1 during add */
      UINTAT(acc+ 4)=0x09090909-UINTAT(acc+ 4);
      UINTAT(acc+ 8)=0x09090909-UINTAT(acc+ 8);
      UINTAT(acc+12)=0x09090909-UINTAT(acc+12);
      UINTAT(acc+16)=0x09090909-UINTAT(acc+16);
      #if QUAD
      UINTAT(acc+20)=0x09090909-UINTAT(acc+20);
      UINTAT(acc+24)=0x09090909-UINTAT(acc+24);
      UINTAT(acc+28)=0x09090909-UINTAT(acc+28);
      UINTAT(acc+32)=0x09090909-UINTAT(acc+32);
      UINTAT(acc+36)=0x09090909-UINTAT(acc+36);
      #endif
      } /* diffsign */
  
    /* now process the rhs coefficient; if it cannot overlap lhs then */
    /* it can be put straight into acc (with an appropriate gap, if */
    /* needed) because no actual addition will be needed (except */
    /* possibly to complete ten's complement) */
    overlap=DECPMAX-(expl-expr);
    #if DECTRACE
    printf("exps: %ld %ld\n", (LI)expl, (LI)expr);
    printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry);
    #endif
  
    if (overlap<=0) {			/* no overlap possible */
      uInt gap;				/* local work */
      /* since a full addition is not needed, a ten's complement */
      /* calculation started above may need to be completed */
      if (carry) {
        for (ub=ulsd; *ub==9; ub--) *ub=0;
        *ub+=1;
        carry=0;				/* taken care of */
        }
      /* up to DECPMAX-1 digits of the final result can extend down */
      /* below the LSD of the lhs, so if the gap is >DECPMAX then the */
      /* rhs will be simply sticky bits.	In this case the gap is */
      /* clamped to DECPMAX and the exponent adjusted to suit [this is */
      /* safe because the lhs is non-zero]. */
      gap=-overlap;
      if (gap>DECPMAX) {
        expr+=gap-1;
        gap=DECPMAX;
        }
      ub=ulsd+gap+1;			/* where MSD will go */
      /* Fill the gap with 0s; note that there is no addition to do */
      ui=&UINTAT(acc+COFF+DECPMAX);	/* start of gap */
      for (; ui<&UINTAT(ub); ui++) *ui=0; /* mind the gap */
      if (overlap<-DECPMAX) {		/* gap was > DECPMAX */
        *ub=(uByte)(!DFISZERO(dfr));	/* make sticky digit */
        }
       else {				/* need full coefficient */
        GETCOEFF(dfr, ub);		/* decode from decFloat */
        ub+=DECPMAX-1;			/* new LSD... */
        }
      ulsd=ub;				/* save new LSD */
      } /* no overlap possible */
  
     else {				/* overlap>0 */
      /* coefficients overlap (perhaps completely, although also */
      /* perhaps only where zeros) */
      ub=buf+COFF+DECPMAX-overlap;	/* where MSD will go */
      /* Fill the prefix gap with 0s; 8 will cover most common */
      /* unalignments, so start with direct assignments (a loop is */
      /* then used for any remaining -- the loop (and the one in a */
      /* moment) is not then on the critical path because the number */
      /* of additions is reduced by (at least) two in this case) */
      UINTAT(buf+4)=0;			/* [clears decQuad 00 too] */
      UINTAT(buf+8)=0;
      if (ub>buf+12) {
        ui=&UINTAT(buf+12);		/* start of any remaining */
        for (; ui<&UINTAT(ub); ui++) *ui=0; /* fill them */
        }
      GETCOEFF(dfr, ub);			/* decode from decFloat */
  
      /* now move tail of rhs across to main acc; again use direct */
      /* assignment for 8 digits-worth */
      UINTAT(acc+COFF+DECPMAX)=UINTAT(buf+COFF+DECPMAX);
      UINTAT(acc+COFF+DECPMAX+4)=UINTAT(buf+COFF+DECPMAX+4);
      if (buf+COFF+DECPMAX+8<ub+DECPMAX) {
        uj=&UINTAT(buf+COFF+DECPMAX+8);	/* source */
        ui=&UINTAT(acc+COFF+DECPMAX+8);	/* target */
        for (; uj<&UINTAT(ub+DECPMAX); ui++, uj++) *ui=*uj;
        }
  
      ulsd=acc+(ub-buf+DECPMAX-1);	/* update LSD pointer */
  
      /* now do the add of the non-tail; this is all nicely aligned, */
      /* and is over a multiple of four digits (because for Quad two */
      /* two 0 digits were added on the left); words in both acc and */
      /* buf (buf especially) will often be zero */
      /* [byte-by-byte add, here, is about 15% slower than the by-fours] */
  
      /* Now effect the add; this is harder on a little-endian */
      /* machine as the inter-digit carry cannot use the usual BCD */
      /* addition trick because the bytes are loaded in the wrong order */
      /* [this loop could be unrolled, but probably scarcely worth it] */
  
      ui=&UINTAT(acc+COFF+DECPMAX-4);	/* target LSW (acc) */
      uj=&UINTAT(buf+COFF+DECPMAX-4);	/* source LSW (buf, to add to acc) */
  
      #if !DECLITEND
      for (; ui>=&UINTAT(acc+4); ui--, uj--) {
        /* bcd8 add */
        carry+=*uj;			/* rhs + carry */
        if (carry==0) continue;		/* no-op */
        carry+=*ui;			/* lhs */
        /* Big-endian BCD adjust (uses internal carry) */
        carry+=0x76f6f6f6;		/* note top nibble not all bits */
        *ui=(carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4); /* BCD adjust */
        carry>>=31;			/* true carry was at far left */
        } /* add loop */
      #else
      for (; ui>=&UINTAT(acc+4); ui--, uj--) {
        /* bcd8 add */
        carry+=*uj;			/* rhs + carry */
        if (carry==0) continue;		/* no-op [common if unaligned] */
        carry+=*ui;			/* lhs */
        /* Little-endian BCD adjust; inter-digit carry must be manual */
        /* because the lsb from the array will be in the most-significant */
        /* byte of carry */
        carry+=0x76767676;		/* note no inter-byte carries */
        carry+=(carry & 0x80000000)>>15;
        carry+=(carry & 0x00800000)>>15;
        carry+=(carry & 0x00008000)>>15;
        carry-=(carry & 0x60606060)>>4;	/* BCD adjust back */
        *ui=carry & 0x0f0f0f0f;		/* clear debris and save */
        /* here, final carry-out bit is at 0x00000080; move it ready */
        /* for next word-add (i.e., to 0x01000000) */
        carry=(carry & 0x00000080)<<17;
        } /* add loop */
      #endif
      #if DECTRACE
      {bcdnum tum;
      printf("Add done, carry=%08lx, diffsign=%ld\n", (LI)carry, (LI)diffsign);
      tum.msd=umsd;  /* acc+4; */
      tum.lsd=ulsd;
      tum.exponent=0;
      tum.sign=0;
      decShowNum(&tum, "dfadd");}
      #endif
      } /* overlap possible */
  
    /* ordering here is a little strange in order to have slowest path */
    /* first in GCC asm listing */
    if (diffsign) {		   /* subtraction */
      if (!carry) {		   /* no carry out means RHS<LHS */
        /* borrowed -- take ten's complement */
        /* sign is lhs sign */
        num.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
  
        /* invert the coefficient first by fours, then add one; space */
        /* at the end of the buffer ensures the by-fours is always */
        /* safe, but lsd+1 must be cleared to prevent a borrow */
        /* if big-endian */
        #if !DECLITEND
        *(ulsd+1)=0;
        #endif
        /* there are always at least four coefficient words */
        UINTAT(umsd)   =0x09090909-UINTAT(umsd);
        UINTAT(umsd+4) =0x09090909-UINTAT(umsd+4);
        UINTAT(umsd+8) =0x09090909-UINTAT(umsd+8);
        UINTAT(umsd+12)=0x09090909-UINTAT(umsd+12);
        #if DOUBLE
  	#define BNEXT 16
        #elif QUAD
  	UINTAT(umsd+16)=0x09090909-UINTAT(umsd+16);
  	UINTAT(umsd+20)=0x09090909-UINTAT(umsd+20);
  	UINTAT(umsd+24)=0x09090909-UINTAT(umsd+24);
  	UINTAT(umsd+28)=0x09090909-UINTAT(umsd+28);
  	UINTAT(umsd+32)=0x09090909-UINTAT(umsd+32);
  	#define BNEXT 36
        #endif
        if (ulsd>=umsd+BNEXT) {		/* unaligned */
  	/* eight will handle most unaligments for Double; 16 for Quad */
  	UINTAT(umsd+BNEXT)=0x09090909-UINTAT(umsd+BNEXT);
  	UINTAT(umsd+BNEXT+4)=0x09090909-UINTAT(umsd+BNEXT+4);
  	#if DOUBLE
  	#define BNEXTY (BNEXT+8)
  	#elif QUAD
  	UINTAT(umsd+BNEXT+8)=0x09090909-UINTAT(umsd+BNEXT+8);
  	UINTAT(umsd+BNEXT+12)=0x09090909-UINTAT(umsd+BNEXT+12);
  	#define BNEXTY (BNEXT+16)
  	#endif
  	if (ulsd>=umsd+BNEXTY) {	/* very unaligned */
  	  ui=&UINTAT(umsd+BNEXTY);	/* -> continue */
  	  for (;;ui++) {
  	    *ui=0x09090909-*ui;		/* invert four digits */
  	    if (ui>=&UINTAT(ulsd-3)) break; /* all done */
  	    }
  	  }
  	}
        /* complete the ten's complement by adding 1 */
        for (ub=ulsd; *ub==9; ub--) *ub=0;
        *ub+=1;
        } /* borrowed */
  
       else {			   /* carry out means RHS>=LHS */
        num.sign=DFWORD(dfr, 0) & DECFLOAT_Sign;
        /* all done except for the special IEEE 754 exact-zero-result */
        /* rule (see above); while testing for zero, strip leading */
        /* zeros (which will save decFinalize doing it) (this is in */
        /* diffsign path, so carry impossible and true umsd is */
        /* acc+COFF) */
  
        /* Check the initial coefficient area using the fast macro; */
        /* this will often be all that needs to be done (as on the */
        /* worst-case path when the subtraction was aligned and */
        /* full-length) */
        if (ISCOEFFZERO(acc+COFF)) {
  	umsd=acc+COFF+DECPMAX-1;   /* so far, so zero */
  	if (ulsd>umsd) {	   /* more to check */
  	  umsd++;		   /* to align after checked area */
  	  for (; UINTAT(umsd)==0 && umsd+3<ulsd;) umsd+=4;
  	  for (; *umsd==0 && umsd<ulsd;) umsd++;
  	  }
  	if (*umsd==0) {		   /* must be true zero (and diffsign) */
  	  num.sign=0;		   /* assume + */
  	  if (set->round==DEC_ROUND_FLOOR) num.sign=DECFLOAT_Sign;
  	  }
  	}
        /* [else was not zero, might still have leading zeros] */
        } /* subtraction gave positive result */
      } /* diffsign */
  
     else { /* same-sign addition */
      num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
      #if DOUBLE
      if (carry) {		   /* only possible with decDouble */
        *(acc+3)=1;		   /* [Quad has leading 00] */
        umsd=acc+3;
        }
      #endif
      } /* same sign */
  
    num.msd=umsd;			   /* set MSD .. */
    num.lsd=ulsd;			   /* .. and LSD */
    num.exponent=expr;		   /* set exponent to smaller */
  
    #if DECTRACE
    decFloatShow(dfl, "dfl");
    decFloatShow(dfr, "dfr");
    decShowNum(&num, "postadd");
    #endif
    return decFinalize(result, &num, set); /* round, check, and lay out */
    } /* decFloatAdd */
  
  /* ------------------------------------------------------------------ */
  /* decFloatAnd -- logical digitwise AND of two decFloats	      */
  /*								      */
  /*   result gets the result of ANDing dfl and dfr		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result, which will be canonical with sign=0	      */
  /*								      */
  /* The operands must be positive, finite with exponent q=0, and	      */
  /* comprise just zeros and ones; if not, Invalid operation results.   */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatAnd(decFloat *result,
  		       const decFloat *dfl, const decFloat *dfr,
  		       decContext *set) {
    if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
     || !DFISCC01(dfl)   || !DFISCC01(dfr)) return decInvalid(result, set);
    /* the operands are positive finite integers (q=0) with just 0s and 1s */
    #if DOUBLE
     DFWORD(result, 0)=ZEROWORD
  		   |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04009124);
     DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x49124491;
    #elif QUAD
     DFWORD(result, 0)=ZEROWORD
  		   |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04000912);
     DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x44912449;
     DFWORD(result, 2)=(DFWORD(dfl, 2) & DFWORD(dfr, 2))&0x12449124;
     DFWORD(result, 3)=(DFWORD(dfl, 3) & DFWORD(dfr, 3))&0x49124491;
    #endif
    return result;
    } /* decFloatAnd */
  
  /* ------------------------------------------------------------------ */
  /* decFloatCanonical -- copy a decFloat, making canonical	      */
  /*								      */
  /*   result gets the canonicalized df				      */
  /*   df	    is the decFloat to copy and make canonical		      */
  /*   returns result						      */
  /*								      */
  /* This works on specials, too; no error or exception is possible.    */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatCanonical(decFloat *result, const decFloat *df) {
    return decCanonical(result, df);
    } /* decFloatCanonical */
  
  /* ------------------------------------------------------------------ */
  /* decFloatClass -- return the class of a decFloat		      */
  /*								      */
  /*   df is the decFloat to test					      */
  /*   returns the decClass that df falls into			      */
  /* ------------------------------------------------------------------ */
  enum decClass decFloatClass(const decFloat *df) {
    Int exp;			   /* exponent */
    if (DFISSPECIAL(df)) {
      if (DFISQNAN(df)) return DEC_CLASS_QNAN;
      if (DFISSNAN(df)) return DEC_CLASS_SNAN;
      /* must be an infinity */
      if (DFISSIGNED(df)) return DEC_CLASS_NEG_INF;
      return DEC_CLASS_POS_INF;
      }
    if (DFISZERO(df)) {		   /* quite common */
      if (DFISSIGNED(df)) return DEC_CLASS_NEG_ZERO;
      return DEC_CLASS_POS_ZERO;
      }
    /* is finite and non-zero; similar code to decFloatIsNormal, here */
    /* [this could be speeded up slightly by in-lining decFloatDigits] */
    exp=GETEXPUN(df)		   /* get unbiased exponent .. */
       +decFloatDigits(df)-1;	   /* .. and make adjusted exponent */
    if (exp>=DECEMIN) {		   /* is normal */
      if (DFISSIGNED(df)) return DEC_CLASS_NEG_NORMAL;
      return DEC_CLASS_POS_NORMAL;
      }
    /* is subnormal */
    if (DFISSIGNED(df)) return DEC_CLASS_NEG_SUBNORMAL;
    return DEC_CLASS_POS_SUBNORMAL;
    } /* decFloatClass */
  
  /* ------------------------------------------------------------------ */
  /* decFloatClassString -- return the class of a decFloat as a string  */
  /*								      */
  /*   df is the decFloat to test					      */
  /*   returns a constant string describing the class df falls into     */
  /* ------------------------------------------------------------------ */
  const char *decFloatClassString(const decFloat *df) {
    enum decClass eclass=decFloatClass(df);
    if (eclass==DEC_CLASS_POS_NORMAL)    return DEC_ClassString_PN;
    if (eclass==DEC_CLASS_NEG_NORMAL)    return DEC_ClassString_NN;
    if (eclass==DEC_CLASS_POS_ZERO)      return DEC_ClassString_PZ;
    if (eclass==DEC_CLASS_NEG_ZERO)      return DEC_ClassString_NZ;
    if (eclass==DEC_CLASS_POS_SUBNORMAL) return DEC_ClassString_PS;
    if (eclass==DEC_CLASS_NEG_SUBNORMAL) return DEC_ClassString_NS;
    if (eclass==DEC_CLASS_POS_INF)       return DEC_ClassString_PI;
    if (eclass==DEC_CLASS_NEG_INF)       return DEC_ClassString_NI;
    if (eclass==DEC_CLASS_QNAN)	       return DEC_ClassString_QN;
    if (eclass==DEC_CLASS_SNAN)	       return DEC_ClassString_SN;
    return DEC_ClassString_UN;	       /* Unknown */
    } /* decFloatClassString */
  
  /* ------------------------------------------------------------------ */
  /* decFloatCompare -- compare two decFloats; quiet NaNs allowed	      */
  /*								      */
  /*   result gets the result of comparing dfl and dfr		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result, which may be -1, 0, 1, or NaN (Unordered)	      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatCompare(decFloat *result,
  			   const decFloat *dfl, const decFloat *dfr,
  			   decContext *set) {
    Int comp;				     /* work */
    /* NaNs are handled as usual */
    if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
    /* numeric comparison needed */
    comp=decNumCompare(dfl, dfr, 0);
    decFloatZero(result);
    if (comp==0) return result;
    DFBYTE(result, DECBYTES-1)=0x01;	/* LSD=1 */
    if (comp<0) DFBYTE(result, 0)|=0x80;	/* set sign bit */
    return result;
    } /* decFloatCompare */
  
  /* ------------------------------------------------------------------ */
  /* decFloatCompareSignal -- compare two decFloats; all NaNs signal    */
  /*								      */
  /*   result gets the result of comparing dfl and dfr		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result, which may be -1, 0, 1, or NaN (Unordered)	      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatCompareSignal(decFloat *result,
  				 const decFloat *dfl, const decFloat *dfr,
  				 decContext *set) {
    Int comp;				     /* work */
    /* NaNs are handled as usual, except that all NaNs signal */
    if (DFISNAN(dfl) || DFISNAN(dfr)) {
      set->status|=DEC_Invalid_operation;
      return decNaNs(result, dfl, dfr, set);
      }
    /* numeric comparison needed */
    comp=decNumCompare(dfl, dfr, 0);
    decFloatZero(result);
    if (comp==0) return result;
    DFBYTE(result, DECBYTES-1)=0x01;	/* LSD=1 */
    if (comp<0) DFBYTE(result, 0)|=0x80;	/* set sign bit */
    return result;
    } /* decFloatCompareSignal */
  
  /* ------------------------------------------------------------------ */
  /* decFloatCompareTotal -- compare two decFloats with total ordering  */
  /*								      */
  /*   result gets the result of comparing dfl and dfr		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   returns result, which may be -1, 0, or 1			      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatCompareTotal(decFloat *result,
  				const decFloat *dfl, const decFloat *dfr) {
    Int comp;				     /* work */
    if (DFISNAN(dfl) || DFISNAN(dfr)) {
      Int nanl, nanr;			     /* work */
      /* morph NaNs to +/- 1 or 2, leave numbers as 0 */
      nanl=DFISSNAN(dfl)+DFISQNAN(dfl)*2;	     /* quiet > signalling */
      if (DFISSIGNED(dfl)) nanl=-nanl;
      nanr=DFISSNAN(dfr)+DFISQNAN(dfr)*2;
      if (DFISSIGNED(dfr)) nanr=-nanr;
      if (nanl>nanr) comp=+1;
       else if (nanl<nanr) comp=-1;
       else { /* NaNs are the same type and sign .. must compare payload */
        /* buffers need +2 for QUAD */
        uByte bufl[DECPMAX+4];		     /* for LHS coefficient + foot */
        uByte bufr[DECPMAX+4];		     /* for RHS coefficient + foot */
        uByte *ub, *uc;			     /* work */
        Int sigl;				     /* signum of LHS */
        sigl=(DFISSIGNED(dfl) ? -1 : +1);
  
        /* decode the coefficients */
        /* (shift both right two if Quad to make a multiple of four) */
        #if QUAD
  	ub = bufl;                           /* avoid type-pun violation */
  	USHORTAT(ub)=0;
  	uc = bufr;                           /* avoid type-pun violation */
  	USHORTAT(uc)=0;
        #endif
        GETCOEFF(dfl, bufl+QUAD*2);	     /* decode from decFloat */
        GETCOEFF(dfr, bufr+QUAD*2);	     /* .. */
        /* all multiples of four, here */
        comp=0;				     /* assume equal */
        for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
  	if (UINTAT(ub)==UINTAT(uc)) continue; /* so far so same */
  	/* about to find a winner; go by bytes in case little-endian */
  	for (;; ub++, uc++) {
  	  if (*ub==*uc) continue;
  	  if (*ub>*uc) comp=sigl;	     /* difference found */
  	   else comp=-sigl;		     /* .. */
  	   break;
  	  }
  	}
        } /* same NaN type and sign */
      }
     else {
      /* numeric comparison needed */
      comp=decNumCompare(dfl, dfr, 1);	/* total ordering */
      }
    decFloatZero(result);
    if (comp==0) return result;
    DFBYTE(result, DECBYTES-1)=0x01;	/* LSD=1 */
    if (comp<0) DFBYTE(result, 0)|=0x80;	/* set sign bit */
    return result;
    } /* decFloatCompareTotal */
  
  /* ------------------------------------------------------------------ */
  /* decFloatCompareTotalMag -- compare magnitudes with total ordering  */
  /*								      */
  /*   result gets the result of comparing abs(dfl) and abs(dfr)	      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   returns result, which may be -1, 0, or 1			      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatCompareTotalMag(decFloat *result,
  				const decFloat *dfl, const decFloat *dfr) {
    decFloat a, b;			/* for copy if needed */
    /* copy and redirect signed operand(s) */
    if (DFISSIGNED(dfl)) {
      decFloatCopyAbs(&a, dfl);
      dfl=&a;
      }
    if (DFISSIGNED(dfr)) {
      decFloatCopyAbs(&b, dfr);
      dfr=&b;
      }
    return decFloatCompareTotal(result, dfl, dfr);
    } /* decFloatCompareTotalMag */
  
  /* ------------------------------------------------------------------ */
  /* decFloatCopy -- copy a decFloat as-is			      */
  /*								      */
  /*   result gets the copy of dfl				      */
  /*   dfl    is the decFloat to copy				      */
  /*   returns result						      */
  /*								      */
  /* This is a bitwise operation; no errors or exceptions are possible. */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) {
    if (dfl!=result) *result=*dfl;	     /* copy needed */
    return result;
    } /* decFloatCopy */
  
  /* ------------------------------------------------------------------ */
  /* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0     */
  /*								      */
  /*   result gets the copy of dfl with sign bit 0		      */
  /*   dfl    is the decFloat to copy				      */
  /*   returns result						      */
  /*								      */
  /* This is a bitwise operation; no errors or exceptions are possible. */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) {
    if (dfl!=result) *result=*dfl;	/* copy needed */
    DFBYTE(result, 0)&=~0x80;		/* zero sign bit */
    return result;
    } /* decFloatCopyAbs */
  
  /* ------------------------------------------------------------------ */
  /* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
  /*								      */
  /*   result gets the copy of dfl with sign bit inverted		      */
  /*   dfl    is the decFloat to copy				      */
  /*   returns result						      */
  /*								      */
  /* This is a bitwise operation; no errors or exceptions are possible. */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) {
    if (dfl!=result) *result=*dfl;	/* copy needed */
    DFBYTE(result, 0)^=0x80;		/* invert sign bit */
    return result;
    } /* decFloatCopyNegate */
  
  /* ------------------------------------------------------------------ */
  /* decFloatCopySign -- copy a decFloat with the sign of another	      */
  /*								      */
  /*   result gets the result of copying dfl with the sign of dfr	      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   returns result						      */
  /*								      */
  /* This is a bitwise operation; no errors or exceptions are possible. */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatCopySign(decFloat *result,
  			    const decFloat *dfl, const decFloat *dfr) {
    uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80);   /* save sign bit */
    if (dfl!=result) *result=*dfl;	     /* copy needed */
    DFBYTE(result, 0)&=~0x80;		     /* clear sign .. */
    DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* .. and set saved */
    return result;
    } /* decFloatCopySign */
  
  /* ------------------------------------------------------------------ */
  /* decFloatDigits -- return the number of digits in a decFloat	      */
  /*								      */
  /*   df is the decFloat to investigate				      */
  /*   returns the number of significant digits in the decFloat; a      */
  /*     zero coefficient returns 1 as does an infinity (a NaN returns  */
  /*     the number of digits in the payload)			      */
  /* ------------------------------------------------------------------ */
  /* private macro to extract a declet according to provided formula */
  /* (form), and if it is non-zero then return the calculated digits */
  /* depending on the declet number (n), where n=0 for the most */
  /* significant declet; uses uInt dpd for work */
  #define dpdlenchk(n, form) {dpd=(form)&0x3ff;	  \
    if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
  /* next one is used when it is known that the declet must be */
  /* non-zero, or is the final zero declet */
  #define dpdlendun(n, form) {dpd=(form)&0x3ff;	  \
    if (dpd==0) return 1;				  \
    return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
  
  uInt decFloatDigits(const decFloat *df) {
    uInt dpd;			   /* work */
    uInt sourhi=DFWORD(df, 0);	   /* top word from source decFloat */
    #if QUAD
    uInt sourmh, sourml;
    #endif
    uInt sourlo;
  
    if (DFISINF(df)) return 1;
    /* A NaN effectively has an MSD of 0; otherwise if non-zero MSD */
    /* then the coefficient is full-length */
    if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX;
  
    #if DOUBLE
      if (sourhi&0x0003ffff) {	 /* ends in first */
        dpdlenchk(0, sourhi>>8);
        sourlo=DFWORD(df, 1);
        dpdlendun(1, (sourhi<<2) | (sourlo>>30));
        } /* [cannot drop through] */
      sourlo=DFWORD(df, 1);  /* sourhi not involved now */
      if (sourlo&0xfff00000) {	 /* in one of first two */
        dpdlenchk(1, sourlo>>30);	 /* very rare */
        dpdlendun(2, sourlo>>20);
        } /* [cannot drop through] */
      dpdlenchk(3, sourlo>>10);
      dpdlendun(4, sourlo);
      /* [cannot drop through] */
  
    #elif QUAD
      if (sourhi&0x00003fff) {	 /* ends in first */
        dpdlenchk(0, sourhi>>4);
        sourmh=DFWORD(df, 1);
        dpdlendun(1, ((sourhi)<<6) | (sourmh>>26));
        } /* [cannot drop through] */
      sourmh=DFWORD(df, 1);
      if (sourmh) {
        dpdlenchk(1, sourmh>>26);
        dpdlenchk(2, sourmh>>16);
        dpdlenchk(3, sourmh>>6);
        sourml=DFWORD(df, 2);
        dpdlendun(4, ((sourmh)<<4) | (sourml>>28));
        } /* [cannot drop through] */
      sourml=DFWORD(df, 2);
      if (sourml) {
        dpdlenchk(4, sourml>>28);
        dpdlenchk(5, sourml>>18);
        dpdlenchk(6, sourml>>8);
        sourlo=DFWORD(df, 3);
        dpdlendun(7, ((sourml)<<2) | (sourlo>>30));
        } /* [cannot drop through] */
      sourlo=DFWORD(df, 3);
      if (sourlo&0xfff00000) {	 /* in one of first two */
        dpdlenchk(7, sourlo>>30);	 /* very rare */
        dpdlendun(8, sourlo>>20);
        } /* [cannot drop through] */
      dpdlenchk(9, sourlo>>10);
      dpdlendun(10, sourlo);
      /* [cannot drop through] */
    #endif
    } /* decFloatDigits */
  
  /* ------------------------------------------------------------------ */
  /* decFloatDivide -- divide a decFloat by another		      */
  /*								      */
  /*   result gets the result of dividing dfl by dfr:		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* ------------------------------------------------------------------ */
  /* This is just a wrapper. */
  decFloat * decFloatDivide(decFloat *result,
  			  const decFloat *dfl, const decFloat *dfr,
  			  decContext *set) {
    return decDivide(result, dfl, dfr, set, DIVIDE);
    } /* decFloatDivide */
  
  /* ------------------------------------------------------------------ */
  /* decFloatDivideInteger -- integer divide a decFloat by another      */
  /*								      */
  /*   result gets the result of dividing dfl by dfr:		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatDivideInteger(decFloat *result,
  			     const decFloat *dfl, const decFloat *dfr,
  			     decContext *set) {
    return decDivide(result, dfl, dfr, set, DIVIDEINT);
    } /* decFloatDivideInteger */
  
  /* ------------------------------------------------------------------ */
  /* decFloatFMA -- multiply and add three decFloats, fused	      */
  /*								      */
  /*   result gets the result of (dfl*dfr)+dff with a single rounding   */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   dff    is the final decFloat (fhs)				      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatFMA(decFloat *result, const decFloat *dfl,
  		       const decFloat *dfr, const decFloat *dff,
  		       decContext *set) {
    /* The accumulator has the bytes needed for FiniteMultiply, plus */
    /* one byte to the left in case of carry, plus DECPMAX+2 to the */
    /* right for the final addition (up to full fhs + round & sticky) */
    #define FMALEN (1+ (DECPMAX9*18) +DECPMAX+2)
    uByte	 acc[FMALEN];		   /* for multiplied coefficient in BCD */
  				   /* .. and for final result */
    bcdnum mul;			   /* for multiplication result */
    bcdnum fin;			   /* for final operand, expanded */
    uByte	 coe[DECPMAX];		   /* dff coefficient in BCD */
    bcdnum *hi, *lo;		   /* bcdnum with higher/lower exponent */
    uInt	 diffsign;		   /* non-zero if signs differ */
    uInt	 hipad;			   /* pad digit for hi if needed */
    Int	 padding;		   /* excess exponent */
    uInt	 carry;			   /* +1 for ten's complement and during add */
    uByte	 *ub, *uh, *ul;		   /* work */
  
    /* handle all the special values [any special operand leads to a */
    /* special result] */
    if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr) || DFISSPECIAL(dff)) {
      decFloat proxy;		   /* multiplication result proxy */
      /* NaNs are handled as usual, giving priority to sNaNs */
      if (DFISSNAN(dfl) || DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
      if (DFISSNAN(dff)) return decNaNs(result, dff, NULL, set);
      if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
      if (DFISNAN(dff)) return decNaNs(result, dff, NULL, set);
      /* One or more of the three is infinite */
      /* infinity times zero is bad */
      decFloatZero(&proxy);
      if (DFISINF(dfl)) {
        if (DFISZERO(dfr)) return decInvalid(result, set);
        decInfinity(&proxy, &proxy);
        }
       else if (DFISINF(dfr)) {
        if (DFISZERO(dfl)) return decInvalid(result, set);
        decInfinity(&proxy, &proxy);
        }
      /* compute sign of multiplication and place in proxy */
      DFWORD(&proxy, 0)|=(DFWORD(dfl, 0)^DFWORD(dfr, 0))&DECFLOAT_Sign;
      if (!DFISINF(dff)) return decFloatCopy(result, &proxy);
      /* dff is Infinite */
      if (!DFISINF(&proxy)) return decInfinity(result, dff);
      /* both sides of addition are infinite; different sign is bad */
      if ((DFWORD(dff, 0)&DECFLOAT_Sign)!=(DFWORD(&proxy, 0)&DECFLOAT_Sign))
        return decInvalid(result, set);
      return decFloatCopy(result, &proxy);
      }
  
    /* Here when all operands are finite */
  
    /* First multiply dfl*dfr */
    decFiniteMultiply(&mul, acc+1, dfl, dfr);
    /* The multiply is complete, exact and unbounded, and described in */
    /* mul with the coefficient held in acc[1...] */
  
    /* now add in dff; the algorithm is essentially the same as */
    /* decFloatAdd, but the code is different because the code there */
    /* is highly optimized for adding two numbers of the same size */
    fin.exponent=GETEXPUN(dff);		/* get dff exponent and sign */
    fin.sign=DFWORD(dff, 0)&DECFLOAT_Sign;
    diffsign=mul.sign^fin.sign;		/* note if signs differ */
    fin.msd=coe;
    fin.lsd=coe+DECPMAX-1;
    GETCOEFF(dff, coe);			/* extract the coefficient */
  
    /* now set hi and lo so that hi points to whichever of mul and fin */
    /* has the higher exponent and lo point to the other [don't care if */
    /* the same] */
    if (mul.exponent>=fin.exponent) {
      hi=&mul;
      lo=&fin;
      }
     else {
      hi=&fin;
      lo=&mul;
      }
  
    /* remove leading zeros on both operands; this will save time later */
    /* and make testing for zero trivial */
    for (; UINTAT(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4;
    for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++;
    for (; UINTAT(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
    for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
  
    /* if hi is zero then result will be lo (which has the smaller */
    /* exponent), which also may need to be tested for zero for the */
    /* weird IEEE 754 sign rules */
    if (*hi->msd==0 && hi->msd==hi->lsd) {     /* hi is zero */
      /* "When the sum of two operands with opposite signs is */
      /* exactly zero, the sign of that sum shall be '+' in all */
      /* rounding modes except round toward -Infinity, in which */
      /* mode that sign shall be '-'." */
      if (diffsign) {
        if (*lo->msd==0 && lo->msd==lo->lsd) { /* lo is zero */
  	lo->sign=0;
  	if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
  	} /* diffsign && lo=0 */
        } /* diffsign */
      return decFinalize(result, lo, set);     /* may need clamping */
      } /* numfl is zero */
    /* [here, both are minimal length and hi is non-zero] */
  
    /* if signs differ, take the ten's complement of hi (zeros to the */
    /* right do not matter because the complement of zero is zero); */
    /* the +1 is done later, as part of the addition, inserted at the */
    /* correct digit */
    hipad=0;
    carry=0;
    if (diffsign) {
      hipad=9;
      carry=1;
      /* exactly the correct number of digits must be inverted */
      for (uh=hi->msd; uh<hi->lsd-3; uh+=4) UINTAT(uh)=0x09090909-UINTAT(uh);
      for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh);
      }
  
    /* ready to add; note that hi has no leading zeros so gap */
    /* calculation does not have to be as pessimistic as in decFloatAdd */
    /* (this is much more like the arbitrary-precision algorithm in */
    /* Rexx and decNumber) */
  
    /* padding is the number of zeros that would need to be added to hi */
    /* for its lsd to be aligned with the lsd of lo */
    padding=hi->exponent-lo->exponent;
    /* printf("FMA pad %ld\n", (LI)padding); */
  
    /* the result of the addition will be built into the accumulator, */
    /* starting from the far right; this could be either hi or lo */
    ub=acc+FMALEN-1;		   /* where lsd of result will go */
    ul=lo->lsd;			   /* lsd of rhs */
  
    if (padding!=0) {		   /* unaligned */
      /* if the msd of lo is more than DECPMAX+2 digits to the right of */
      /* the original msd of hi then it can be reduced to a single */
      /* digit at the right place, as it stays clear of hi digits */
      /* [it must be DECPMAX+2 because during a subtraction the msd */
      /* could become 0 after a borrow from 1.000 to 0.9999...] */
      Int hilen=(Int)(hi->lsd-hi->msd+1); /* lengths */
      Int lolen=(Int)(lo->lsd-lo->msd+1); /* .. */
      Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
      Int reduce=newexp-lo->exponent;
      if (reduce>0) {			/* [= case gives reduce=0 nop] */
        /* printf("FMA reduce: %ld\n", (LI)reduce); */
        if (reduce>=lolen) {		/* eating all */
  	lo->lsd=lo->msd;		/* reduce to single digit */
  	lo->exponent=newexp;		/* [known to be non-zero] */
  	}
         else { /* < */
  	uByte *up=lo->lsd;
  	lo->lsd=lo->lsd-reduce;
  	if (*lo->lsd==0)		/* could need sticky bit */
  	 for (; up>lo->lsd; up--) {	/* search discarded digits */
  	  if (*up!=0) {			/* found one... */
  	    *lo->lsd=1;			/* set sticky bit */
  	    break;
  	    }
  	  }
  	lo->exponent+=reduce;
  	}
        padding=hi->exponent-lo->exponent; /* recalculate */
        ul=lo->lsd;			 /* .. */
        } /* maybe reduce */
      /* padding is now <= DECPMAX+2 but still > 0; tricky DOUBLE case */
      /* is when hi is a 1 that will become a 0.9999... by subtraction: */
      /*	 hi:   1				   E+16 */
      /*	 lo:	.................1000000000000000  E-16 */
      /* which for the addition pads and reduces to: */
      /*	 hi:   1000000000000000000		   E-2 */
      /*	 lo:	.................1		   E-2 */
      #if DECCHECK
        if (padding>DECPMAX+2) printf("FMA excess padding: %ld\n", (LI)padding);
        if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding);
        /* printf("FMA padding: %ld\n", (LI)padding); */
      #endif
      /* padding digits can now be set in the result; one or more of */
      /* these will come from lo; others will be zeros in the gap */
      for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul;
      for (;padding>0; padding--, ub--) *ub=0; /* mind the gap */
      }
  
    /* addition now complete to the right of the rightmost digit of hi */
    uh=hi->lsd;
  
    /* carry was set up depending on ten's complement above; do the add... */
    for (;; ub--) {
      uInt hid, lod;
      if (uh<hi->msd) {
        if (ul<lo->msd) break;
        hid=hipad;
        }
       else hid=*uh--;
      if (ul<lo->msd) lod=0;
       else lod=*ul--;
      *ub=(uByte)(carry+hid+lod);
      if (*ub<10) carry=0;
       else {
        *ub-=10;
        carry=1;
        }
      } /* addition loop */
  
    /* addition complete -- now handle carry, borrow, etc. */
    /* use lo to set up the num (its exponent is already correct, and */
    /* sign usually is) */
    lo->msd=ub+1;
    lo->lsd=acc+FMALEN-1;
    /* decShowNum(lo, "lo"); */
    if (!diffsign) {		   /* same-sign addition */
      if (carry) {		   /* carry out */
        *ub=1;			   /* place the 1 .. */
        lo->msd--;		   /* .. and update */
        }
      } /* same sign */
     else {			   /* signs differed (subtraction) */
      if (!carry) {		   /* no carry out means hi<lo */
        /* borrowed -- take ten's complement of the right digits */
        lo->sign=hi->sign;	   /* sign is lhs sign */
        for (ul=lo->msd; ul<lo->lsd-3; ul+=4) UINTAT(ul)=0x09090909-UINTAT(ul);
        for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); /* [leaves ul at lsd+1] */
        /* complete the ten's complement by adding 1 [cannot overrun] */
        for (ul--; *ul==9; ul--) *ul=0;
        *ul+=1;
        } /* borrowed */
       else {			   /* carry out means hi>=lo */
        /* sign to use is lo->sign */
        /* all done except for the special IEEE 754 exact-zero-result */
        /* rule (see above); while testing for zero, strip leading */
        /* zeros (which will save decFinalize doing it) */
        for (; UINTAT(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
        for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
        if (*lo->msd==0) {	   /* must be true zero (and diffsign) */
  	lo->sign=0;		   /* assume + */
  	if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
  	}
        /* [else was not zero, might still have leading zeros] */
        } /* subtraction gave positive result */
      } /* diffsign */
  
    return decFinalize(result, lo, set);	/* round, check, and lay out */
    } /* decFloatFMA */
  
  /* ------------------------------------------------------------------ */
  /* decFloatFromInt -- initialise a decFloat from an Int		      */
  /*								      */
  /*   result gets the converted Int				      */
  /*   n	    is the Int to convert				      */
  /*   returns result						      */
  /*								      */
  /* The result is Exact; no errors or exceptions are possible.	      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatFromInt32(decFloat *result, Int n) {
    uInt u=(uInt)n;			/* copy as bits */
    uInt encode;				/* work */
    DFWORD(result, 0)=ZEROWORD;		/* always */
    #if QUAD
      DFWORD(result, 1)=0;
      DFWORD(result, 2)=0;
    #endif
    if (n<0) {				/* handle -n with care */
      /* [This can be done without the test, but is then slightly slower] */
      u=(~u)+1;
      DFWORD(result, 0)|=DECFLOAT_Sign;
      }
    /* Since the maximum value of u now is 2**31, only the low word of */
    /* result is affected */
    encode=BIN2DPD[u%1000];
    u/=1000;
    encode|=BIN2DPD[u%1000]<<10;
    u/=1000;
    encode|=BIN2DPD[u%1000]<<20;
    u/=1000;				/* now 0, 1, or 2 */
    encode|=u<<30;
    DFWORD(result, DECWORDS-1)=encode;
    return result;
    } /* decFloatFromInt32 */
  
  /* ------------------------------------------------------------------ */
  /* decFloatFromUInt -- initialise a decFloat from a uInt	      */
  /*								      */
  /*   result gets the converted uInt				      */
  /*   n	    is the uInt to convert				      */
  /*   returns result						      */
  /*								      */
  /* The result is Exact; no errors or exceptions are possible.	      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatFromUInt32(decFloat *result, uInt u) {
    uInt encode;				/* work */
    DFWORD(result, 0)=ZEROWORD;		/* always */
    #if QUAD
      DFWORD(result, 1)=0;
      DFWORD(result, 2)=0;
    #endif
    encode=BIN2DPD[u%1000];
    u/=1000;
    encode|=BIN2DPD[u%1000]<<10;
    u/=1000;
    encode|=BIN2DPD[u%1000]<<20;
    u/=1000;				/* now 0 -> 4 */
    encode|=u<<30;
    DFWORD(result, DECWORDS-1)=encode;
    DFWORD(result, DECWORDS-2)|=u>>2;	/* rarely non-zero */
    return result;
    } /* decFloatFromUInt32 */
  
  /* ------------------------------------------------------------------ */
  /* decFloatInvert -- logical digitwise INVERT of a decFloat	      */
  /*								      */
  /*   result gets the result of INVERTing df			      */
  /*   df	    is the decFloat to invert				      */
  /*   set    is the context					      */
  /*   returns result, which will be canonical with sign=0	      */
  /*								      */
  /* The operand must be positive, finite with exponent q=0, and	      */
  /* comprise just zeros and ones; if not, Invalid operation results.   */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatInvert(decFloat *result, const decFloat *df,
  			  decContext *set) {
    uInt sourhi=DFWORD(df, 0);		/* top word of dfs */
  
    if (!DFISUINT01(df) || !DFISCC01(df)) return decInvalid(result, set);
    /* the operand is a finite integer (q=0) */
    #if DOUBLE
     DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04009124);
     DFWORD(result, 1)=(~DFWORD(df, 1))	&0x49124491;
    #elif QUAD
     DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04000912);
     DFWORD(result, 1)=(~DFWORD(df, 1))	&0x44912449;
     DFWORD(result, 2)=(~DFWORD(df, 2))	&0x12449124;
     DFWORD(result, 3)=(~DFWORD(df, 3))	&0x49124491;
    #endif
    return result;
    } /* decFloatInvert */
  
  /* ------------------------------------------------------------------ */
  /* decFloatIs -- decFloat tests (IsSigned, etc.)		      */
  /*								      */
  /*   df is the decFloat to test					      */
  /*   returns 0 or 1 in an int32_t				      */
  /*								      */
  /* Many of these could be macros, but having them as real functions   */
  /* is a bit cleaner (and they can be referred to here by the generic  */
  /* names)							      */
  /* ------------------------------------------------------------------ */
  uInt decFloatIsCanonical(const decFloat *df) {
    if (DFISSPECIAL(df)) {
      if (DFISINF(df)) {
        if (DFWORD(df, 0)&ECONMASK) return 0;  /* exponent continuation */
        if (!DFISCCZERO(df)) return 0;	     /* coefficient continuation */
        return 1;
        }
      /* is a NaN */
      if (DFWORD(df, 0)&ECONNANMASK) return 0; /* exponent continuation */
      if (DFISCCZERO(df)) return 1;	     /* coefficient continuation */
      /* drop through to check payload */
      }
    { /* declare block */
    #if DOUBLE
      uInt sourhi=DFWORD(df, 0);
      uInt sourlo=DFWORD(df, 1);
      if (CANONDPDOFF(sourhi, 8)
       && CANONDPDTWO(sourhi, sourlo, 30)
       && CANONDPDOFF(sourlo, 20)
       && CANONDPDOFF(sourlo, 10)
       && CANONDPDOFF(sourlo, 0)) return 1;
    #elif QUAD
      uInt sourhi=DFWORD(df, 0);
      uInt sourmh=DFWORD(df, 1);
      uInt sourml=DFWORD(df, 2);
      uInt sourlo=DFWORD(df, 3);
      if (CANONDPDOFF(sourhi, 4)
       && CANONDPDTWO(sourhi, sourmh, 26)
       && CANONDPDOFF(sourmh, 16)
       && CANONDPDOFF(sourmh, 6)
       && CANONDPDTWO(sourmh, sourml, 28)
       && CANONDPDOFF(sourml, 18)
       && CANONDPDOFF(sourml, 8)
       && CANONDPDTWO(sourml, sourlo, 30)
       && CANONDPDOFF(sourlo, 20)
       && CANONDPDOFF(sourlo, 10)
       && CANONDPDOFF(sourlo, 0)) return 1;
    #endif
    } /* block */
    return 0;    /* a declet is non-canonical */
    }
  
  uInt decFloatIsFinite(const decFloat *df) {
    return !DFISSPECIAL(df);
    }
  uInt decFloatIsInfinite(const decFloat *df) {
    return DFISINF(df);
    }
  uInt decFloatIsInteger(const decFloat *df) {
    return DFISINT(df);
    }
  uInt decFloatIsNaN(const decFloat *df) {
    return DFISNAN(df);
    }
  uInt decFloatIsNormal(const decFloat *df) {
    Int exp;			   /* exponent */
    if (DFISSPECIAL(df)) return 0;
    if (DFISZERO(df)) return 0;
    /* is finite and non-zero */
    exp=GETEXPUN(df)		   /* get unbiased exponent .. */
       +decFloatDigits(df)-1;	   /* .. and make adjusted exponent */
    return (exp>=DECEMIN);	   /* < DECEMIN is subnormal */
    }
  uInt decFloatIsSignaling(const decFloat *df) {
    return DFISSNAN(df);
    }
  uInt decFloatIsSignalling(const decFloat *df) {
    return DFISSNAN(df);
    }
  uInt decFloatIsSigned(const decFloat *df) {
    return DFISSIGNED(df);
    }
  uInt decFloatIsSubnormal(const decFloat *df) {
    if (DFISSPECIAL(df)) return 0;
    /* is finite */
    if (decFloatIsNormal(df)) return 0;
    /* it is <Nmin, but could be zero */
    if (DFISZERO(df)) return 0;
    return 1;				     /* is subnormal */
    }
  uInt decFloatIsZero(const decFloat *df) {
    return DFISZERO(df);
    } /* decFloatIs... */
  
  /* ------------------------------------------------------------------ */
  /* decFloatLogB -- return adjusted exponent, by 754r rules	      */
  /*								      */
  /*   result gets the adjusted exponent as an integer, or a NaN etc.   */
  /*   df	    is the decFloat to be examined			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* Notable cases:						      */
  /*   A<0 -> Use |A|						      */
  /*   A=0 -> -Infinity (Division by zero)			      */
  /*   A=Infinite -> +Infinity (Exact)				      */
  /*   A=1 exactly -> 0 (Exact)					      */
  /*   NaNs are propagated as usual				      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatLogB(decFloat *result, const decFloat *df,
  			decContext *set) {
    Int ae;				     /* adjusted exponent */
    if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
    if (DFISINF(df)) {
      DFWORD(result, 0)=0;		     /* need +ve */
      return decInfinity(result, result);	     /* canonical +Infinity */
      }
    if (DFISZERO(df)) {
      set->status|=DEC_Division_by_zero;	     /* as per 754r */
      DFWORD(result, 0)=DECFLOAT_Sign;	     /* make negative */
      return decInfinity(result, result);	     /* canonical -Infinity */
      }
    ae=GETEXPUN(df)			/* get unbiased exponent .. */
      +decFloatDigits(df)-1;		/* .. and make adjusted exponent */
    /* ae has limited range (3 digits for DOUBLE and 4 for QUAD), so */
    /* it is worth using a special case of decFloatFromInt32 */
    DFWORD(result, 0)=ZEROWORD;		/* always */
    if (ae<0) {
      DFWORD(result, 0)|=DECFLOAT_Sign;	/* -0 so far */
      ae=-ae;
      }
    #if DOUBLE
      DFWORD(result, 1)=BIN2DPD[ae];	/* a single declet */
    #elif QUAD
      DFWORD(result, 1)=0;
      DFWORD(result, 2)=0;
      DFWORD(result, 3)=(ae/1000)<<10;	/* is <10, so need no DPD encode */
      DFWORD(result, 3)|=BIN2DPD[ae%1000];
    #endif
    return result;
    } /* decFloatLogB */
  
  /* ------------------------------------------------------------------ */
  /* decFloatMax -- return maxnum of two operands			      */
  /*								      */
  /*   result gets the chosen decFloat				      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* If just one operand is a quiet NaN it is ignored.		      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatMax(decFloat *result,
  		       const decFloat *dfl, const decFloat *dfr,
  		       decContext *set) {
    Int comp;
    if (DFISNAN(dfl)) {
      /* sNaN or both NaNs leads to normal NaN processing */
      if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
      return decCanonical(result, dfr);	     /* RHS is numeric */
      }
    if (DFISNAN(dfr)) {
      /* sNaN leads to normal NaN processing (both NaN handled above) */
      if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
      return decCanonical(result, dfl);	     /* LHS is numeric */
      }
    /* Both operands are numeric; numeric comparison needed -- use */
    /* total order for a well-defined choice (and +0 > -0) */
    comp=decNumCompare(dfl, dfr, 1);
    if (comp>=0) return decCanonical(result, dfl);
    return decCanonical(result, dfr);
    } /* decFloatMax */
  
  /* ------------------------------------------------------------------ */
  /* decFloatMaxMag -- return maxnummag of two operands		      */
  /*								      */
  /*   result gets the chosen decFloat				      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* Returns according to the magnitude comparisons if both numeric and */
  /* unequal, otherwise returns maxnum				      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatMaxMag(decFloat *result,
  		       const decFloat *dfl, const decFloat *dfr,
  		       decContext *set) {
    Int comp;
    decFloat absl, absr;
    if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMax(result, dfl, dfr, set);
  
    decFloatCopyAbs(&absl, dfl);
    decFloatCopyAbs(&absr, dfr);
    comp=decNumCompare(&absl, &absr, 0);
    if (comp>0) return decCanonical(result, dfl);
    if (comp<0) return decCanonical(result, dfr);
    return decFloatMax(result, dfl, dfr, set);
    } /* decFloatMaxMag */
  
  /* ------------------------------------------------------------------ */
  /* decFloatMin -- return minnum of two operands			      */
  /*								      */
  /*   result gets the chosen decFloat				      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* If just one operand is a quiet NaN it is ignored.		      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatMin(decFloat *result,
  		       const decFloat *dfl, const decFloat *dfr,
  		       decContext *set) {
    Int comp;
    if (DFISNAN(dfl)) {
      /* sNaN or both NaNs leads to normal NaN processing */
      if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
      return decCanonical(result, dfr);	     /* RHS is numeric */
      }
    if (DFISNAN(dfr)) {
      /* sNaN leads to normal NaN processing (both NaN handled above) */
      if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
      return decCanonical(result, dfl);	     /* LHS is numeric */
      }
    /* Both operands are numeric; numeric comparison needed -- use */
    /* total order for a well-defined choice (and +0 > -0) */
    comp=decNumCompare(dfl, dfr, 1);
    if (comp<=0) return decCanonical(result, dfl);
    return decCanonical(result, dfr);
    } /* decFloatMin */
  
  /* ------------------------------------------------------------------ */
  /* decFloatMinMag -- return minnummag of two operands		      */
  /*								      */
  /*   result gets the chosen decFloat				      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* Returns according to the magnitude comparisons if both numeric and */
  /* unequal, otherwise returns minnum				      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatMinMag(decFloat *result,
  		       const decFloat *dfl, const decFloat *dfr,
  		       decContext *set) {
    Int comp;
    decFloat absl, absr;
    if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMin(result, dfl, dfr, set);
  
    decFloatCopyAbs(&absl, dfl);
    decFloatCopyAbs(&absr, dfr);
    comp=decNumCompare(&absl, &absr, 0);
    if (comp<0) return decCanonical(result, dfl);
    if (comp>0) return decCanonical(result, dfr);
    return decFloatMin(result, dfl, dfr, set);
    } /* decFloatMinMag */
  
  /* ------------------------------------------------------------------ */
  /* decFloatMinus -- negate value, heeding NaNs, etc.		      */
  /*								      */
  /*   result gets the canonicalized 0-df				      */
  /*   df	    is the decFloat to minus				      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* This has the same effect as 0-df where the exponent of the zero is */
  /* the same as that of df (if df is finite).			      */
  /* The effect is also the same as decFloatCopyNegate except that NaNs */
  /* are handled normally (the sign of a NaN is not affected, and an    */
  /* sNaN will signal), the result is canonical, and zero gets sign 0.  */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatMinus(decFloat *result, const decFloat *df,
  			 decContext *set) {
    if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
    decCanonical(result, df);			  /* copy and check */
    if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80;	  /* turn off sign bit */
     else DFBYTE(result, 0)^=0x80;		  /* flip sign bit */
    return result;
    } /* decFloatMinus */
  
  /* ------------------------------------------------------------------ */
  /* decFloatMultiply -- multiply two decFloats			      */
  /*								      */
  /*   result gets the result of multiplying dfl and dfr:		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatMultiply(decFloat *result,
  			    const decFloat *dfl, const decFloat *dfr,
  			    decContext *set) {
    bcdnum num;			   /* for final conversion */
    uByte	 bcdacc[DECPMAX9*18+1];	   /* for coefficent in BCD */
  
    if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
      /* NaNs are handled as usual */
      if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
      /* infinity times zero is bad */
      if (DFISINF(dfl) && DFISZERO(dfr)) return decInvalid(result, set);
      if (DFISINF(dfr) && DFISZERO(dfl)) return decInvalid(result, set);
      /* both infinite; return canonical infinity with computed sign */
      DFWORD(result, 0)=DFWORD(dfl, 0)^DFWORD(dfr, 0); /* compute sign */
      return decInfinity(result, result);
      }
  
    /* Here when both operands are finite */
    decFiniteMultiply(&num, bcdacc, dfl, dfr);
    return decFinalize(result, &num, set); /* round, check, and lay out */
    } /* decFloatMultiply */
  
  /* ------------------------------------------------------------------ */
  /* decFloatNextMinus -- next towards -Infinity			      */
  /*								      */
  /*   result gets the next lesser decFloat			      */
  /*   dfl    is the decFloat to start with			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* This is 754r nextdown; Invalid is the only status possible (from   */
  /* an sNaN).							      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
  			     decContext *set) {
    decFloat delta;			/* tiny increment */
    uInt savestat;			/* saves status */
    enum rounding saveround;		/* .. and mode */
  
    /* +Infinity is the special case */
    if (DFISINF(dfl) && !DFISSIGNED(dfl)) {
      DFSETNMAX(result);
      return result;			/* [no status to set] */
      }
    /* other cases are effected by sutracting a tiny delta -- this */
    /* should be done in a wider format as the delta is unrepresentable */
    /* here (but can be done with normal add if the sign of zero is */
    /* treated carefully, because no Inexactitude is interesting); */
    /* rounding to -Infinity then pushes the result to next below */
    decFloatZero(&delta);			/* set up tiny delta */
    DFWORD(&delta, DECWORDS-1)=1;		/* coefficient=1 */
    DFWORD(&delta, 0)=DECFLOAT_Sign;	/* Sign=1 + biased exponent=0 */
    /* set up for the directional round */
    saveround=set->round;			/* save mode */
    set->round=DEC_ROUND_FLOOR;		/* .. round towards -Infinity */
    savestat=set->status;			/* save status */
    decFloatAdd(result, dfl, &delta, set);
    /* Add rules mess up the sign when going from +Ntiny to 0 */
    if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
    set->status&=DEC_Invalid_operation;	/* preserve only sNaN status */
    set->status|=savestat;		/* restore pending flags */
    set->round=saveround;			/* .. and mode */
    return result;
    } /* decFloatNextMinus */
  
  /* ------------------------------------------------------------------ */
  /* decFloatNextPlus -- next towards +Infinity			      */
  /*								      */
  /*   result gets the next larger decFloat			      */
  /*   dfl    is the decFloat to start with			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* This is 754r nextup; Invalid is the only status possible (from     */
  /* an sNaN).							      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
  			    decContext *set) {
    uInt savestat;			/* saves status */
    enum rounding saveround;		/* .. and mode */
    decFloat delta;			/* tiny increment */
  
    /* -Infinity is the special case */
    if (DFISINF(dfl) && DFISSIGNED(dfl)) {
      DFSETNMAX(result);
      DFWORD(result, 0)|=DECFLOAT_Sign;	/* make negative */
      return result;			/* [no status to set] */
      }
    /* other cases are effected by sutracting a tiny delta -- this */
    /* should be done in a wider format as the delta is unrepresentable */
    /* here (but can be done with normal add if the sign of zero is */
    /* treated carefully, because no Inexactitude is interesting); */
    /* rounding to +Infinity then pushes the result to next above */
    decFloatZero(&delta);			/* set up tiny delta */
    DFWORD(&delta, DECWORDS-1)=1;		/* coefficient=1 */
    DFWORD(&delta, 0)=0;			/* Sign=0 + biased exponent=0 */
    /* set up for the directional round */
    saveround=set->round;			/* save mode */
    set->round=DEC_ROUND_CEILING;		/* .. round towards +Infinity */
    savestat=set->status;			/* save status */
    decFloatAdd(result, dfl, &delta, set);
    /* Add rules mess up the sign when going from -Ntiny to -0 */
    if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
    set->status&=DEC_Invalid_operation;	/* preserve only sNaN status */
    set->status|=savestat;		/* restore pending flags */
    set->round=saveround;			/* .. and mode */
    return result;
    } /* decFloatNextPlus */
  
  /* ------------------------------------------------------------------ */
  /* decFloatNextToward -- next towards a decFloat		      */
  /*								      */
  /*   result gets the next decFloat				      */
  /*   dfl    is the decFloat to start with			      */
  /*   dfr    is the decFloat to move toward			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* This is 754r nextafter; status may be set unless the result is a   */
  /* normal number.						      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatNextToward(decFloat *result,
  			      const decFloat *dfl, const decFloat *dfr,
  			      decContext *set) {
    decFloat delta;			/* tiny increment or decrement */
    decFloat pointone;			/* 1e-1 */
    uInt	savestat;			/* saves status */
    enum	rounding saveround;		/* .. and mode */
    uInt	deltatop;			/* top word for delta */
    Int	comp;				/* work */
  
    if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
    /* Both are numeric, so Invalid no longer a possibility */
    comp=decNumCompare(dfl, dfr, 0);
    if (comp==0) return decFloatCopySign(result, dfl, dfr); /* equal */
    /* unequal; do NextPlus or NextMinus but with different status rules */
  
    if (comp<0) { /* lhs<rhs, do NextPlus, see above for commentary */
      if (DFISINF(dfl) && DFISSIGNED(dfl)) {   /* -Infinity special case */
        DFSETNMAX(result);
        DFWORD(result, 0)|=DECFLOAT_Sign;
        return result;
        }
      saveround=set->round;		     /* save mode */
      set->round=DEC_ROUND_CEILING;	     /* .. round towards +Infinity */
      deltatop=0;				     /* positive delta */
      }
     else { /* lhs>rhs, do NextMinus, see above for commentary */
      if (DFISINF(dfl) && !DFISSIGNED(dfl)) {  /* +Infinity special case */
        DFSETNMAX(result);
        return result;
        }
      saveround=set->round;		     /* save mode */
      set->round=DEC_ROUND_FLOOR;		     /* .. round towards -Infinity */
      deltatop=DECFLOAT_Sign;		     /* negative delta */
      }
    savestat=set->status;			     /* save status */
    /* Here, Inexact is needed where appropriate (and hence Underflow, */
    /* etc.).  Therefore the tiny delta which is otherwise */
    /* unrepresentable (see NextPlus and NextMinus) is constructed */
    /* using the multiplication of FMA. */
    decFloatZero(&delta);			/* set up tiny delta */
    DFWORD(&delta, DECWORDS-1)=1;		/* coefficient=1 */
    DFWORD(&delta, 0)=deltatop;		/* Sign + biased exponent=0 */
    decFloatFromString(&pointone, "1E-1", set); /* set up multiplier */
    decFloatFMA(result, &delta, &pointone, dfl, set);
    /* [Delta is truly tiny, so no need to correct sign of zero] */
    /* use new status unless the result is normal */
    if (decFloatIsNormal(result)) set->status=savestat; /* else goes forward */
    set->round=saveround;			/* restore mode */
    return result;
    } /* decFloatNextToward */
  
  /* ------------------------------------------------------------------ */
  /* decFloatOr -- logical digitwise OR of two decFloats		      */
  /*								      */
  /*   result gets the result of ORing dfl and dfr		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result, which will be canonical with sign=0	      */
  /*								      */
  /* The operands must be positive, finite with exponent q=0, and	      */
  /* comprise just zeros and ones; if not, Invalid operation results.   */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatOr(decFloat *result,
  		       const decFloat *dfl, const decFloat *dfr,
  		       decContext *set) {
    if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
     || !DFISCC01(dfl)   || !DFISCC01(dfr)) return decInvalid(result, set);
    /* the operands are positive finite integers (q=0) with just 0s and 1s */
    #if DOUBLE
     DFWORD(result, 0)=ZEROWORD
  		   |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04009124);
     DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x49124491;
    #elif QUAD
     DFWORD(result, 0)=ZEROWORD
  		   |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04000912);
     DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x44912449;
     DFWORD(result, 2)=(DFWORD(dfl, 2) | DFWORD(dfr, 2))&0x12449124;
     DFWORD(result, 3)=(DFWORD(dfl, 3) | DFWORD(dfr, 3))&0x49124491;
    #endif
    return result;
    } /* decFloatOr */
  
  /* ------------------------------------------------------------------ */
  /* decFloatPlus -- add value to 0, heeding NaNs, etc.		      */
  /*								      */
  /*   result gets the canonicalized 0+df				      */
  /*   df	    is the decFloat to plus				      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* This has the same effect as 0+df where the exponent of the zero is */
  /* the same as that of df (if df is finite).			      */
  /* The effect is also the same as decFloatCopy except that NaNs	      */
  /* are handled normally (the sign of a NaN is not affected, and an    */
  /* sNaN will signal), the result is canonical, and zero gets sign 0.  */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatPlus(decFloat *result, const decFloat *df,
  			decContext *set) {
    if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
    decCanonical(result, df);			  /* copy and check */
    if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80;	  /* turn off sign bit */
    return result;
    } /* decFloatPlus */
  
  /* ------------------------------------------------------------------ */
  /* decFloatQuantize -- quantize a decFloat			      */
  /*								      */
  /*   result gets the result of quantizing dfl to match dfr	      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs), which sets the exponent     */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* Unless there is an error or the result is infinite, the exponent   */
  /* of result is guaranteed to be the same as that of dfr.	      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatQuantize(decFloat *result,
  			    const decFloat *dfl, const decFloat *dfr,
  			    decContext *set) {
    Int	explb, exprb;	      /* left and right biased exponents */
    uByte *ulsd;		      /* local LSD pointer */
    uInt	*ui;		      /* work */
    uByte *ub;		      /* .. */
    Int	drop;		      /* .. */
    uInt	dpd;		      /* .. */
    uInt	encode;		      /* encoding accumulator */
    uInt	sourhil, sourhir;     /* top words from source decFloats */
    /* the following buffer holds the coefficient for manipulation */
    uByte buf[4+DECPMAX*3];     /* + space for zeros to left or right */
    #if DECTRACE
    bcdnum num;		      /* for trace displays */
    #endif
  
    /* Start decoding the arguments */
    sourhil=DFWORD(dfl, 0);	   /* LHS top word */
    explb=DECCOMBEXP[sourhil>>26];   /* get exponent high bits (in place) */
    sourhir=DFWORD(dfr, 0);	   /* RHS top word */
    exprb=DECCOMBEXP[sourhir>>26];
  
    if (EXPISSPECIAL(explb | exprb)) { /* either is special? */
      /* NaNs are handled as usual */
      if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
      /* one infinity but not both is bad */
      if (DFISINF(dfl)!=DFISINF(dfr)) return decInvalid(result, set);
      /* both infinite; return canonical infinity with sign of LHS */
      return decInfinity(result, dfl);
      }
  
    /* Here when both arguments are finite */
    /* complete extraction of the exponents [no need to unbias] */
    explb+=GETECON(dfl);		   /* + continuation */
    exprb+=GETECON(dfr);		   /* .. */
  
    /* calculate the number of digits to drop from the coefficient */
    drop=exprb-explb;		   /* 0 if nothing to do */
    if (drop==0) return decCanonical(result, dfl); /* return canonical */
  
    /* the coefficient is needed; lay it out into buf, offset so zeros */
    /* can be added before or after as needed -- an extra heading is */
    /* added so can safely pad Quad DECPMAX-1 zeros to the left by */
    /* fours */
    #define BUFOFF (buf+4+DECPMAX)
    GETCOEFF(dfl, BUFOFF);	   /* decode from decFloat */
    /* [now the msd is at BUFOFF and the lsd is at BUFOFF+DECPMAX-1] */
  
    #if DECTRACE
    num.msd=BUFOFF;
    num.lsd=BUFOFF+DECPMAX-1;
    num.exponent=explb-DECBIAS;
    num.sign=sourhil & DECFLOAT_Sign;
    decShowNum(&num, "dfl");
    #endif
  
    if (drop>0) {				/* [most common case] */
      /* (this code is very similar to that in decFloatFinalize, but */
      /* has many differences so is duplicated here -- so any changes */
      /* may need to be made there, too) */
      uByte *roundat;			     /* -> re-round digit */
      uByte reround;			     /* reround value */
      /* printf("Rounding; drop=%ld\n", (LI)drop); */
  
      /* there is at least one zero needed to the left, in all but one */
      /* exceptional (all-nines) case, so place four zeros now; this is */
      /* needed almost always and makes rounding all-nines by fours safe */
      UINTAT(BUFOFF-4)=0;
  
      /* Three cases here: */
      /*	 1. new LSD is in coefficient (almost always) */
      /*	 2. new LSD is digit to left of coefficient (so MSD is */
      /*	    round-for-reround digit) */
      /*	 3. new LSD is to left of case 2 (whole coefficient is sticky) */
      /* Note that leading zeros can safely be treated as useful digits */
  
      /* [duplicate check-stickies code to save a test] */
      /* [by-digit check for stickies as runs of zeros are rare] */
      if (drop<DECPMAX) {			     /* NB lengths not addresses */
        roundat=BUFOFF+DECPMAX-drop;
        reround=*roundat;
        for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
  	if (*ub!=0) {			     /* non-zero to be discarded */
  	  reround=DECSTICKYTAB[reround];     /* apply sticky bit */
  	  break;			     /* [remainder don't-care] */
  	  }
  	} /* check stickies */
        ulsd=roundat-1;			     /* set LSD */
        }
       else {				     /* edge case */
        if (drop==DECPMAX) {
  	roundat=BUFOFF;
  	reround=*roundat;
  	}
         else {
  	roundat=BUFOFF-1;
  	reround=0;
  	}
        for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
  	if (*ub!=0) {			     /* non-zero to be discarded */
  	  reround=DECSTICKYTAB[reround];     /* apply sticky bit */
  	  break;			     /* [remainder don't-care] */
  	  }
  	} /* check stickies */
        *BUFOFF=0;			     /* make a coefficient of 0 */
        ulsd=BUFOFF;			     /* .. at the MSD place */
        }
  
      if (reround!=0) {			     /* discarding non-zero */
        uInt bump=0;
        set->status|=DEC_Inexact;
  
        /* next decide whether to increment the coefficient */
        if (set->round==DEC_ROUND_HALF_EVEN) { /* fastpath slowest case */
  	if (reround>5) bump=1;		     /* >0.5 goes up */
  	 else if (reround==5)		     /* exactly 0.5000 .. */
  	  bump=*ulsd & 0x01;		     /* .. up iff [new] lsd is odd */
  	} /* r-h-e */
         else switch (set->round) {
  	case DEC_ROUND_DOWN: {
  	  /* no change */
  	  break;} /* r-d */
  	case DEC_ROUND_HALF_DOWN: {
  	  if (reround>5) bump=1;
  	  break;} /* r-h-d */
  	case DEC_ROUND_HALF_UP: {
  	  if (reround>=5) bump=1;
  	  break;} /* r-h-u */
  	case DEC_ROUND_UP: {
  	  if (reround>0) bump=1;
  	  break;} /* r-u */
  	case DEC_ROUND_CEILING: {
  	  /* same as _UP for positive numbers, and as _DOWN for negatives */
  	  if (!(sourhil&DECFLOAT_Sign) && reround>0) bump=1;
  	  break;} /* r-c */
  	case DEC_ROUND_FLOOR: {
  	  /* same as _UP for negative numbers, and as _DOWN for positive */
  	  /* [negative reround cannot occur on 0] */
  	  if (sourhil&DECFLOAT_Sign && reround>0) bump=1;
  	  break;} /* r-f */
  	case DEC_ROUND_05UP: {
  	  if (reround>0) { /* anything out there is 'sticky' */
  	    /* bump iff lsd=0 or 5; this cannot carry so it could be */
  	    /* effected immediately with no bump -- but the code */
  	    /* is clearer if this is done the same way as the others */
  	    if (*ulsd==0 || *ulsd==5) bump=1;
  	    }
  	  break;} /* r-r */
  	default: {	/* e.g., DEC_ROUND_MAX */
  	  set->status|=DEC_Invalid_context;
  	  #if DECCHECK
  	  printf("Unknown rounding mode: %ld\n", (LI)set->round);
  	  #endif
  	  break;}
  	} /* switch (not r-h-e) */
        /* printf("ReRound: %ld  bump: %ld\n", (LI)reround, (LI)bump); */
  
        if (bump!=0) {			     /* need increment */
  	/* increment the coefficient; this could give 1000... (after */
  	/* the all nines case) */
  	ub=ulsd;
  	for (; UINTAT(ub-3)==0x09090909; ub-=4) UINTAT(ub-3)=0;
  	/* now at most 3 digits left to non-9 (usually just the one) */
  	for (; *ub==9; ub--) *ub=0;
  	*ub+=1;
  	/* [the all-nines case will have carried one digit to the */
  	/* left of the original MSD -- just where it is needed] */
  	} /* bump needed */
        } /* inexact rounding */
  
      /* now clear zeros to the left so exactly DECPMAX digits will be */
      /* available in the coefficent -- the first word to the left was */
      /* cleared earlier for safe carry; now add any more needed */
      if (drop>4) {
        UINTAT(BUFOFF-8)=0;		     /* must be at least 5 */
        for (ui=&UINTAT(BUFOFF-12); ui>&UINTAT(ulsd-DECPMAX-3); ui--) *ui=0;
        }
      } /* need round (drop>0) */
  
     else { /* drop<0; padding with -drop digits is needed */
      /* This is the case where an error can occur if the padded */
      /* coefficient will not fit; checking for this can be done in the */
      /* same loop as padding for zeros if the no-hope and zero cases */
      /* are checked first */
      if (-drop>DECPMAX-1) {		     /* cannot fit unless 0 */
        if (!ISCOEFFZERO(BUFOFF)) return decInvalid(result, set);
        /* a zero can have any exponent; just drop through and use it */
        ulsd=BUFOFF+DECPMAX-1;
        }
       else { /* padding will fit (but may still be too long) */
        /* final-word mask depends on endianess */
        #if DECLITEND
        static const uInt dmask[]={0, 0x000000ff, 0x0000ffff, 0x00ffffff};
        #else
        static const uInt dmask[]={0, 0xff000000, 0xffff0000, 0xffffff00};
        #endif
        for (ui=&UINTAT(BUFOFF+DECPMAX);; ui++) {
  	*ui=0;
  	if (UINTAT(&UBYTEAT(ui)-DECPMAX)!=0) { /* could be bad */
  	  /* if all four digits should be zero, definitely bad */
  	  if (ui<=&UINTAT(BUFOFF+DECPMAX+(-drop)-4))
  	    return decInvalid(result, set);
  	  /* must be a 1- to 3-digit sequence; check more carefully */
  	  if ((UINTAT(&UBYTEAT(ui)-DECPMAX)&dmask[(-drop)%4])!=0)
  	    return decInvalid(result, set);
  	  break;    /* no need for loop end test */
  	  }
  	if (ui>=&UINTAT(BUFOFF+DECPMAX+(-drop)-4)) break; /* done */
  	}
        ulsd=BUFOFF+DECPMAX+(-drop)-1;
        } /* pad and check leading zeros */
      } /* drop<0 */
  
    #if DECTRACE
    num.msd=ulsd-DECPMAX+1;
    num.lsd=ulsd;
    num.exponent=explb-DECBIAS;
    num.sign=sourhil & DECFLOAT_Sign;
    decShowNum(&num, "res");
    #endif
  
    /*------------------------------------------------------------------*/
    /* At this point the result is DECPMAX digits, ending at ulsd, so   */
    /* fits the encoding exactly; there is no possibility of error      */
    /*------------------------------------------------------------------*/
    encode=((exprb>>DECECONL)<<4) + *(ulsd-DECPMAX+1); /* make index */
    encode=DECCOMBFROM[encode];		     /* indexed by (0-2)*16+msd */
    /* the exponent continuation can be extracted from the original RHS */
    encode|=sourhir & ECONMASK;
    encode|=sourhil&DECFLOAT_Sign;	     /* add the sign from LHS */
  
    /* finally encode the coefficient */
    /* private macro to encode a declet; this version can be used */
    /* because all coefficient digits exist */
    #define getDPD3q(dpd, n) ub=ulsd-(3*(n))-2;			\
      dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];
  
    #if DOUBLE
      getDPD3q(dpd, 4); encode|=dpd<<8;
      getDPD3q(dpd, 3); encode|=dpd>>2;
      DFWORD(result, 0)=encode;
      encode=dpd<<30;
      getDPD3q(dpd, 2); encode|=dpd<<20;
      getDPD3q(dpd, 1); encode|=dpd<<10;
      getDPD3q(dpd, 0); encode|=dpd;
      DFWORD(result, 1)=encode;
  
    #elif QUAD
      getDPD3q(dpd,10); encode|=dpd<<4;
      getDPD3q(dpd, 9); encode|=dpd>>6;
      DFWORD(result, 0)=encode;
      encode=dpd<<26;
      getDPD3q(dpd, 8); encode|=dpd<<16;
      getDPD3q(dpd, 7); encode|=dpd<<6;
      getDPD3q(dpd, 6); encode|=dpd>>4;
      DFWORD(result, 1)=encode;
      encode=dpd<<28;
      getDPD3q(dpd, 5); encode|=dpd<<18;
      getDPD3q(dpd, 4); encode|=dpd<<8;
      getDPD3q(dpd, 3); encode|=dpd>>2;
      DFWORD(result, 2)=encode;
      encode=dpd<<30;
      getDPD3q(dpd, 2); encode|=dpd<<20;
      getDPD3q(dpd, 1); encode|=dpd<<10;
      getDPD3q(dpd, 0); encode|=dpd;
      DFWORD(result, 3)=encode;
    #endif
    return result;
    } /* decFloatQuantize */
  
  /* ------------------------------------------------------------------ */
  /* decFloatReduce -- reduce finite coefficient to minimum length      */
  /*								      */
  /*   result gets the reduced decFloat				      */
  /*   df	    is the source decFloat				      */
  /*   set    is the context					      */
  /*   returns result, which will be canonical			      */
  /*								      */
  /* This removes all possible trailing zeros from the coefficient;     */
  /* some may remain when the number is very close to Nmax.	      */
  /* Special values are unchanged and no status is set unless df=sNaN.  */
  /* Reduced zero has an exponent q=0.				      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatReduce(decFloat *result, const decFloat *df,
  			  decContext *set) {
    bcdnum num;				/* work */
    uByte buf[DECPMAX], *ub;		/* coefficient and pointer */
    if (df!=result) *result=*df;		/* copy, if needed */
    if (DFISNAN(df)) return decNaNs(result, df, NULL, set);   /* sNaN */
    /* zeros and infinites propagate too */
    if (DFISINF(df)) return decInfinity(result, df);     /* canonical */
    if (DFISZERO(df)) {
      uInt sign=DFWORD(df, 0)&DECFLOAT_Sign;
      decFloatZero(result);
      DFWORD(result, 0)|=sign;
      return result;			/* exponent dropped, sign OK */
      }
    /* non-zero finite */
    GETCOEFF(df, buf);
    ub=buf+DECPMAX-1;			/* -> lsd */
    if (*ub) return result;		/* no trailing zeros */
    for (ub--; *ub==0;) ub--;		/* terminates because non-zero */
    /* *ub is the first non-zero from the right */
    num.sign=DFWORD(df, 0)&DECFLOAT_Sign; /* set up number... */
    num.exponent=GETEXPUN(df)+(Int)(buf+DECPMAX-1-ub); /* adjusted exponent */
    num.msd=buf;
    num.lsd=ub;
    return decFinalize(result, &num, set);
    } /* decFloatReduce */
  
  /* ------------------------------------------------------------------ */
  /* decFloatRemainder -- integer divide and return remainder	      */
  /*								      */
  /*   result gets the remainder of dividing dfl by dfr:		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatRemainder(decFloat *result,
  			     const decFloat *dfl, const decFloat *dfr,
  			     decContext *set) {
    return decDivide(result, dfl, dfr, set, REMAINDER);
    } /* decFloatRemainder */
  
  /* ------------------------------------------------------------------ */
  /* decFloatRemainderNear -- integer divide to nearest and remainder   */
  /*								      */
  /*   result gets the remainder of dividing dfl by dfr:		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* This is the IEEE remainder, where the nearest integer is used.     */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatRemainderNear(decFloat *result,
  			     const decFloat *dfl, const decFloat *dfr,
  			     decContext *set) {
    return decDivide(result, dfl, dfr, set, REMNEAR);
    } /* decFloatRemainderNear */
  
  /* ------------------------------------------------------------------ */
  /* decFloatRotate -- rotate the coefficient of a decFloat left/right  */
  /*								      */
  /*   result gets the result of rotating dfl			      */
  /*   dfl    is the source decFloat to rotate			      */
  /*   dfr    is the count of digits to rotate, an integer (with q=0)   */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* The digits of the coefficient of dfl are rotated to the left (if   */
  /* dfr is positive) or to the right (if dfr is negative) without      */
  /* adjusting the exponent or the sign of dfl.			      */
  /*								      */
  /* dfr must be in the range -DECPMAX through +DECPMAX.		      */
  /* NaNs are propagated as usual.  An infinite dfl is unaffected (but  */
  /* dfr must be valid).	No status is set unless dfr is invalid or an  */
  /* operand is an sNaN.	The result is canonical.		      */
  /* ------------------------------------------------------------------ */
  #define PHALF (ROUNDUP(DECPMAX/2, 4))	/* half length, rounded up */
  decFloat * decFloatRotate(decFloat *result,
  			 const decFloat *dfl, const decFloat *dfr,
  			 decContext *set) {
    Int rotate;				/* dfr as an Int */
    uByte buf[DECPMAX+PHALF];		/* coefficient + half */
    uInt digits, savestat;		/* work */
    bcdnum num;				/* .. */
    uByte *ub;				/* .. */
  
    if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
    if (!DFISINT(dfr)) return decInvalid(result, set);
    digits=decFloatDigits(dfr);			 /* calculate digits */
    if (digits>2) return decInvalid(result, set);	 /* definitely out of range */
    rotate=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
    if (rotate>DECPMAX) return decInvalid(result, set); /* too big */
    /* [from here on no error or status change is possible] */
    if (DFISINF(dfl)) return decInfinity(result, dfl);  /* canonical */
    /* handle no-rotate cases */
    if (rotate==0 || rotate==DECPMAX) return decCanonical(result, dfl);
    /* a real rotate is needed: 0 < rotate < DECPMAX */
    /* reduce the rotation to no more than half to reduce copying later */
    /* (for QUAD in fact half + 2 digits) */
    if (DFISSIGNED(dfr)) rotate=-rotate;
    if (abs(rotate)>PHALF) {
      if (rotate<0) rotate=DECPMAX+rotate;
       else rotate=rotate-DECPMAX;
      }
    /* now lay out the coefficient, leaving room to the right or the */
    /* left depending on the direction of rotation */
    ub=buf;
    if (rotate<0) ub+=PHALF;    /* rotate right, so space to left */
    GETCOEFF(dfl, ub);
    /* copy half the digits to left or right, and set num.msd */
    if (rotate<0) {
      memcpy(buf, buf+DECPMAX, PHALF);
      num.msd=buf+PHALF+rotate;
      }
     else {
      memcpy(buf+DECPMAX, buf, PHALF);
      num.msd=buf+rotate;
      }
    /* fill in rest of num */
    num.lsd=num.msd+DECPMAX-1;
    num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
    num.exponent=GETEXPUN(dfl);
    savestat=set->status;			/* record */
    decFinalize(result, &num, set);
    set->status=savestat;			/* restore */
    return result;
    } /* decFloatRotate */
  
  /* ------------------------------------------------------------------ */
  /* decFloatSameQuantum -- test decFloats for same quantum	      */
  /*								      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   returns 1 if the operands have the same quantum, 0 otherwise     */
  /*								      */
  /* No error is possible and no status results.			      */
  /* ------------------------------------------------------------------ */
  uInt decFloatSameQuantum(const decFloat *dfl, const decFloat *dfr) {
    if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) {
      if (DFISNAN(dfl) && DFISNAN(dfr)) return 1;
      if (DFISINF(dfl) && DFISINF(dfr)) return 1;
      return 0;  /* any other special mixture gives false */
      }
    if (GETEXP(dfl)==GETEXP(dfr)) return 1; /* biased exponents match */
    return 0;
    } /* decFloatSameQuantum */
  
  /* ------------------------------------------------------------------ */
  /* decFloatScaleB -- multiply by a power of 10, as per 754r	      */
  /*								      */
  /*   result gets the result of the operation			      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs), am integer (with q=0)	      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* This computes result=dfl x 10**dfr where dfr is an integer in the  */
  /* range +/-2*(emax+pmax), typically resulting from LogB.	      */
  /* Underflow and Overflow (with Inexact) may occur.  NaNs propagate   */
  /* as usual.							      */
  /* ------------------------------------------------------------------ */
  #define SCALEBMAX 2*(DECEMAX+DECPMAX)	/* D=800, Q=12356 */
  decFloat * decFloatScaleB(decFloat *result,
  			  const decFloat *dfl, const decFloat *dfr,
  			  decContext *set) {
    uInt digits;				/* work */
    Int  expr;				/* dfr as an Int */
  
    if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
    if (!DFISINT(dfr)) return decInvalid(result, set);
    digits=decFloatDigits(dfr);		     /* calculate digits */
  
    #if DOUBLE
    if (digits>3) return decInvalid(result, set);	  /* definitely out of range */
    expr=DPD2BIN[DFWORD(dfr, 1)&0x3ff];		  /* must be in bottom declet */
    #elif QUAD
    if (digits>5) return decInvalid(result, set);	  /* definitely out of range */
    expr=DPD2BIN[DFWORD(dfr, 3)&0x3ff]		  /* in bottom 2 declets .. */
        +DPD2BIN[(DFWORD(dfr, 3)>>10)&0x3ff]*1000;  /* .. */
    #endif
    if (expr>SCALEBMAX) return decInvalid(result, set);  /* oops */
    /* [from now on no error possible] */
    if (DFISINF(dfl)) return decInfinity(result, dfl);   /* canonical */
    if (DFISSIGNED(dfr)) expr=-expr;
    /* dfl is finite and expr is valid */
    *result=*dfl;				     /* copy to target */
    return decFloatSetExponent(result, set, GETEXPUN(result)+expr);
    } /* decFloatScaleB */
  
  /* ------------------------------------------------------------------ */
  /* decFloatShift -- shift the coefficient of a decFloat left or right */
  /*								      */
  /*   result gets the result of shifting dfl			      */
  /*   dfl    is the source decFloat to shift			      */
  /*   dfr    is the count of digits to shift, an integer (with q=0)    */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* The digits of the coefficient of dfl are shifted to the left (if   */
  /* dfr is positive) or to the right (if dfr is negative) without      */
  /* adjusting the exponent or the sign of dfl.			      */
  /*								      */
  /* dfr must be in the range -DECPMAX through +DECPMAX.		      */
  /* NaNs are propagated as usual.  An infinite dfl is unaffected (but  */
  /* dfr must be valid).	No status is set unless dfr is invalid or an  */
  /* operand is an sNaN.	The result is canonical.		      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatShift(decFloat *result,
  			 const decFloat *dfl, const decFloat *dfr,
  			 decContext *set) {
    Int shift;				/* dfr as an Int */
    uByte buf[DECPMAX*2];			/* coefficient + padding */
    uInt digits, savestat;		/* work */
    bcdnum num;				/* .. */
  
    if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
    if (!DFISINT(dfr)) return decInvalid(result, set);
    digits=decFloatDigits(dfr);			  /* calculate digits */
    if (digits>2) return decInvalid(result, set);	  /* definitely out of range */
    shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff];	  /* is in bottom declet */
    if (shift>DECPMAX) return decInvalid(result, set);   /* too big */
    /* [from here on no error or status change is possible] */
  
    if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
    /* handle no-shift and all-shift (clear to zero) cases */
    if (shift==0) return decCanonical(result, dfl);
    if (shift==DECPMAX) {			     /* zero with sign */
      uByte sign=(uByte)(DFBYTE(dfl, 0)&0x80); /* save sign bit */
      decFloatZero(result);		     /* make +0 */
      DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* and set sign */
      /* [cannot safely use CopySign] */
      return result;
      }
    /* a real shift is needed: 0 < shift < DECPMAX */
    num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
    num.exponent=GETEXPUN(dfl);
    num.msd=buf;
    GETCOEFF(dfl, buf);
    if (DFISSIGNED(dfr)) { /* shift right */
      /* edge cases are taken care of, so this is easy */
      num.lsd=buf+DECPMAX-shift-1;
      }
     else { /* shift left -- zero padding needed to right */
      UINTAT(buf+DECPMAX)=0;		/* 8 will handle most cases */
      UINTAT(buf+DECPMAX+4)=0;		/* .. */
      if (shift>8) memset(buf+DECPMAX+8, 0, 8+QUAD*18); /* all other cases */
      num.msd+=shift;
      num.lsd=num.msd+DECPMAX-1;
      }
    savestat=set->status;			/* record */
    decFinalize(result, &num, set);
    set->status=savestat;			/* restore */
    return result;
    } /* decFloatShift */
  
  /* ------------------------------------------------------------------ */
  /* decFloatSubtract -- subtract a decFloat from another		      */
  /*								      */
  /*   result gets the result of subtracting dfr from dfl:	      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatSubtract(decFloat *result,
  			    const decFloat *dfl, const decFloat *dfr,
  			    decContext *set) {
    decFloat temp;
    /* NaNs must propagate without sign change */
    if (DFISNAN(dfr)) return decFloatAdd(result, dfl, dfr, set);
    temp=*dfr;				       /* make a copy */
    DFBYTE(&temp, 0)^=0x80;		       /* flip sign */
    return decFloatAdd(result, dfl, &temp, set); /* and add to the lhs */
    } /* decFloatSubtract */
  
  /* ------------------------------------------------------------------ */
  /* decFloatToInt -- round to 32-bit binary integer (4 flavours)	      */
  /*								      */
  /*   df	   is the decFloat to round				      */
  /*   set   is the context					      */
  /*   round is the rounding mode to use				      */
  /*   returns a uInt or an Int, rounded according to the name	      */
  /*								      */
  /* Invalid will always be signaled if df is a NaN, is Infinite, or is */
  /* outside the range of the target; Inexact will not be signaled for  */
  /* simple rounding unless 'Exact' appears in the name.		      */
  /* ------------------------------------------------------------------ */
  uInt decFloatToUInt32(const decFloat *df, decContext *set,
  		      enum rounding round) {
    return decToInt32(df, set, round, 0, 1);}
  
  uInt decFloatToUInt32Exact(const decFloat *df, decContext *set,
  			   enum rounding round) {
    return decToInt32(df, set, round, 1, 1);}
  
  Int decFloatToInt32(const decFloat *df, decContext *set,
  		    enum rounding round) {
    return (Int)decToInt32(df, set, round, 0, 0);}
  
  Int decFloatToInt32Exact(const decFloat *df, decContext *set,
  			 enum rounding round) {
    return (Int)decToInt32(df, set, round, 1, 0);}
  
  /* ------------------------------------------------------------------ */
  /* decFloatToIntegral -- round to integral value (two flavours)	      */
  /*								      */
  /*   result gets the result					      */
  /*   df	    is the decFloat to round				      */
  /*   set    is the context					      */
  /*   round  is the rounding mode to use				      */
  /*   returns result						      */
  /*								      */
  /* No exceptions, even Inexact, are raised except for sNaN input, or  */
  /* if 'Exact' appears in the name.				      */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatToIntegralValue(decFloat *result, const decFloat *df,
  				   decContext *set, enum rounding round) {
    return decToIntegral(result, df, set, round, 0);}
  
  decFloat * decFloatToIntegralExact(decFloat *result, const decFloat *df,
  				   decContext *set) {
    return decToIntegral(result, df, set, set->round, 1);}
  
  /* ------------------------------------------------------------------ */
  /* decFloatXor -- logical digitwise XOR of two decFloats	      */
  /*								      */
  /*   result gets the result of XORing dfl and dfr		      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs)			      */
  /*   set    is the context					      */
  /*   returns result, which will be canonical with sign=0	      */
  /*								      */
  /* The operands must be positive, finite with exponent q=0, and	      */
  /* comprise just zeros and ones; if not, Invalid operation results.   */
  /* ------------------------------------------------------------------ */
  decFloat * decFloatXor(decFloat *result,
  		       const decFloat *dfl, const decFloat *dfr,
  		       decContext *set) {
    if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
     || !DFISCC01(dfl)   || !DFISCC01(dfr)) return decInvalid(result, set);
    /* the operands are positive finite integers (q=0) with just 0s and 1s */
    #if DOUBLE
     DFWORD(result, 0)=ZEROWORD
  		   |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04009124);
     DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x49124491;
    #elif QUAD
     DFWORD(result, 0)=ZEROWORD
  		   |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04000912);
     DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x44912449;
     DFWORD(result, 2)=(DFWORD(dfl, 2) ^ DFWORD(dfr, 2))&0x12449124;
     DFWORD(result, 3)=(DFWORD(dfl, 3) ^ DFWORD(dfr, 3))&0x49124491;
    #endif
    return result;
    } /* decFloatXor */
  
  /* ------------------------------------------------------------------ */
  /* decInvalid -- set Invalid_operation result			      */
  /*								      */
  /*   result gets a canonical NaN				      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* status has Invalid_operation added				      */
  /* ------------------------------------------------------------------ */
  static decFloat *decInvalid(decFloat *result, decContext *set) {
    decFloatZero(result);
    DFWORD(result, 0)=DECFLOAT_qNaN;
    set->status|=DEC_Invalid_operation;
    return result;
    } /* decInvalid */
  
  /* ------------------------------------------------------------------ */
  /* decInfinity -- set canonical Infinity with sign from a decFloat    */
  /*								      */
  /*   result gets a canonical Infinity				      */
  /*   df	    is source decFloat (only the sign is used)		      */
  /*   returns result						      */
  /*								      */
  /* df may be the same as result					      */
  /* ------------------------------------------------------------------ */
  static decFloat *decInfinity(decFloat *result, const decFloat *df) {
    uInt sign=DFWORD(df, 0);	   /* save source signword */
    decFloatZero(result);		   /* clear everything */
    DFWORD(result, 0)=DECFLOAT_Inf | (sign & DECFLOAT_Sign);
    return result;
    } /* decInfinity */
  
  /* ------------------------------------------------------------------ */
  /* decNaNs -- handle NaN argument(s)				      */
  /*								      */
  /*   result gets the result of handling dfl and dfr, one or both of   */
  /*	    which is a NaN					      */
  /*   dfl    is the first decFloat (lhs)				      */
  /*   dfr    is the second decFloat (rhs) -- may be NULL for a single- */
  /*	    operand operation					      */
  /*   set    is the context					      */
  /*   returns result						      */
  /*								      */
  /* Called when one or both operands is a NaN, and propagates the      */
  /* appropriate result to res.  When an sNaN is found, it is changed   */
  /* to a qNaN and Invalid operation is set.			      */
  /* ------------------------------------------------------------------ */
  static decFloat *decNaNs(decFloat *result,
  			 const decFloat *dfl, const decFloat *dfr,
  			 decContext *set) {
    /* handle sNaNs first */
    if (dfr!=NULL && DFISSNAN(dfr) && !DFISSNAN(dfl)) dfl=dfr; /* use RHS */
    if (DFISSNAN(dfl)) {
      decCanonical(result, dfl);		/* propagate canonical sNaN */
      DFWORD(result, 0)&=~(DECFLOAT_qNaN ^ DECFLOAT_sNaN); /* quiet */
      set->status|=DEC_Invalid_operation;
      return result;
      }
    /* one or both is a quiet NaN */
    if (!DFISNAN(dfl)) dfl=dfr;		/* RHS must be NaN, use it */
    return decCanonical(result, dfl);	/* propagate canonical qNaN */
    } /* decNaNs */
  
  /* ------------------------------------------------------------------ */
  /* decNumCompare -- numeric comparison of two decFloats		      */
  /*								      */
  /*   dfl    is the left-hand decFloat, which is not a NaN	      */
  /*   dfr    is the right-hand decFloat, which is not a NaN	      */
  /*   tot    is 1 for total order compare, 0 for simple numeric	      */
  /*   returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr		      */
  /*								      */
  /* No error is possible; status and mode are unchanged.		      */
  /* ------------------------------------------------------------------ */
  static Int decNumCompare(const decFloat *dfl, const decFloat *dfr, Flag tot) {
    Int	sigl, sigr;			/* LHS and RHS non-0 signums */
    Int	shift;				/* shift needed to align operands */
    uByte *ub, *uc;			/* work */
    /* buffers +2 if Quad (36 digits), need double plus 4 for safe padding */
    uByte bufl[DECPMAX*2+QUAD*2+4];	/* for LHS coefficient + padding */
    uByte bufr[DECPMAX*2+QUAD*2+4];	/* for RHS coefficient + padding */
  
    sigl=1;
    if (DFISSIGNED(dfl)) {
      if (!DFISSIGNED(dfr)) {		/* -LHS +RHS */
        if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
        return -1;			/* RHS wins */
        }
      sigl=-1;
      }
    if (DFISSIGNED(dfr)) {
      if (!DFISSIGNED(dfl)) {		/* +LHS -RHS */
        if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
        return +1;			/* LHS wins */
        }
      }
  
    /* signs are the same; operand(s) could be zero */
    sigr=-sigl;				/* sign to return if abs(RHS) wins */
  
    if (DFISINF(dfl)) {
      if (DFISINF(dfr)) return 0;		/* both infinite & same sign */
      return sigl;			/* inf > n */
      }
    if (DFISINF(dfr)) return sigr;	/* n < inf [dfl is finite] */
  
    /* here, both are same sign and finite; calculate their offset */
    shift=GETEXP(dfl)-GETEXP(dfr);	/* [0 means aligned] */
    /* [bias can be ignored -- the absolute exponent is not relevant] */
  
    if (DFISZERO(dfl)) {
      if (!DFISZERO(dfr)) return sigr;	/* LHS=0, RHS!=0 */
      /* both are zero, return 0 if both same exponent or numeric compare */
      if (shift==0 || !tot) return 0;
      if (shift>0) return sigl;
      return sigr;			/* [shift<0] */
      }
     else {				/* LHS!=0 */
      if (DFISZERO(dfr)) return sigl;	/* LHS!=0, RHS=0 */
      }
    /* both are known to be non-zero at this point */
  
    /* if the exponents are so different that the coefficients do not */
    /* overlap (by even one digit) then a full comparison is not needed */
    if (abs(shift)>=DECPMAX) {		/* no overlap */
      /* coefficients are known to be non-zero */
      if (shift>0) return sigl;
      return sigr;			/* [shift<0] */
      }
  
    /* decode the coefficients */
    /* (shift both right two if Quad to make a multiple of four) */
    #if QUAD
      ub=bufl;                            /* avoid type-pun violation */
      UINTAT(ub)=0;
      uc=bufr;                            /* avoid type-pun violation */
      UINTAT(uc)=0;
    #endif
    GETCOEFF(dfl, bufl+QUAD*2);		/* decode from decFloat */
    GETCOEFF(dfr, bufr+QUAD*2);		/* .. */
    if (shift==0) {			/* aligned; common and easy */
      /* all multiples of four, here */
      for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
        if (UINTAT(ub)==UINTAT(uc)) continue; /* so far so same */
        /* about to find a winner; go by bytes in case little-endian */
        for (;; ub++, uc++) {
  	if (*ub>*uc) return sigl;	/* difference found */
  	if (*ub<*uc) return sigr;	/* .. */
  	}
        }
      } /* aligned */
     else if (shift>0) {			/* lhs to left */
      ub=bufl;				/* RHS pointer */
      /* pad bufl so right-aligned; most shifts will fit in 8 */
      UINTAT(bufl+DECPMAX+QUAD*2)=0;	/* add eight zeros */
      UINTAT(bufl+DECPMAX+QUAD*2+4)=0;	/* .. */
      if (shift>8) {
        /* more than eight; fill the rest, and also worth doing the */
        /* lead-in by fours */
        uByte *up;			 /* work */
        uByte *upend=bufl+DECPMAX+QUAD*2+shift;
        for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
        /* [pads up to 36 in all for Quad] */
        for (;; ub+=4) {
  	if (UINTAT(ub)!=0) return sigl;
  	if (ub+4>bufl+shift-4) break;
  	}
        }
      /* check remaining leading digits */
      for (; ub<bufl+shift; ub++) if (*ub!=0) return sigl;
      /* now start the overlapped part; bufl has been padded, so the */
      /* comparison can go for the full length of bufr, which is a */
      /* multiple of 4 bytes */
      for (uc=bufr; ; uc+=4, ub+=4) {
        if (UINTAT(uc)!=UINTAT(ub)) {	/* mismatch found */
  	for (;; uc++, ub++) {		/* check from left [little-endian?] */
  	  if (*ub>*uc) return sigl;	/* difference found */
  	  if (*ub<*uc) return sigr;	/* .. */
  	  }
  	} /* mismatch */
        if (uc==bufr+QUAD*2+DECPMAX-4) break; /* all checked */
        }
      } /* shift>0 */
  
     else { /* shift<0) .. RHS is to left of LHS; mirror shift>0 */
      uc=bufr;				/* RHS pointer */
      /* pad bufr so right-aligned; most shifts will fit in 8 */
      UINTAT(bufr+DECPMAX+QUAD*2)=0;	/* add eight zeros */
      UINTAT(bufr+DECPMAX+QUAD*2+4)=0;	/* .. */
      if (shift<-8) {
        /* more than eight; fill the rest, and also worth doing the */
        /* lead-in by fours */
        uByte *up;			 /* work */
        uByte *upend=bufr+DECPMAX+QUAD*2-shift;
        for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
        /* [pads up to 36 in all for Quad] */
        for (;; uc+=4) {
  	if (UINTAT(uc)!=0) return sigr;
  	if (uc+4>bufr-shift-4) break;
  	}
        }
      /* check remaining leading digits */
      for (; uc<bufr-shift; uc++) if (*uc!=0) return sigr;
      /* now start the overlapped part; bufr has been padded, so the */
      /* comparison can go for the full length of bufl, which is a */
      /* multiple of 4 bytes */
      for (ub=bufl; ; ub+=4, uc+=4) {
        if (UINTAT(ub)!=UINTAT(uc)) {	/* mismatch found */
  	for (;; ub++, uc++) {		/* check from left [little-endian?] */
  	  if (*ub>*uc) return sigl;	/* difference found */
  	  if (*ub<*uc) return sigr;	/* .. */
  	  }
  	} /* mismatch */
        if (ub==bufl+QUAD*2+DECPMAX-4) break; /* all checked */
        }
      } /* shift<0 */
  
    /* Here when compare equal */
    if (!tot) return 0;			/* numerically equal */
    /* total ordering .. exponent matters */
    if (shift>0) return sigl;		/* total order by exponent */
    if (shift<0) return sigr;		/* .. */
    return 0;
    } /* decNumCompare */
  
  /* ------------------------------------------------------------------ */
  /* decToInt32 -- local routine to effect ToInteger conversions	      */
  /*								      */
  /*   df	    is the decFloat to convert				      */
  /*   set    is the context					      */
  /*   rmode  is the rounding mode to use				      */
  /*   exact  is 1 if Inexact should be signalled			      */
  /*   unsign is 1 if the result a uInt, 0 if an Int (cast to uInt)     */
  /*   returns 32-bit result as a uInt				      */
  /*								      */
  /* Invalid is set is df is a NaN, is infinite, or is out-of-range; in */
  /* these cases 0 is returned.					      */
  /* ------------------------------------------------------------------ */
  static uInt decToInt32(const decFloat *df, decContext *set,
  		       enum rounding rmode, Flag exact, Flag unsign) {
    Int  exp;			   /* exponent */
    uInt sourhi, sourpen, sourlo;	   /* top word from source decFloat .. */
    uInt hi, lo;			   /* .. penultimate, least, etc. */
    decFloat zero, result;	   /* work */
    Int  i;			   /* .. */
  
    /* Start decoding the argument */
    sourhi=DFWORD(df, 0);			/* top word */
    exp=DECCOMBEXP[sourhi>>26];		/* get exponent high bits (in place) */
    if (EXPISSPECIAL(exp)) {		/* is special? */
      set->status|=DEC_Invalid_operation; /* signal */
      return 0;
      }
  
    /* Here when the argument is finite */
    if (GETEXPUN(df)==0) result=*df;	/* already a true integer */
     else {				/* need to round to integer */
      enum rounding saveround;		/* saver */
      uInt savestatus;			/* .. */
      saveround=set->round;		/* save rounding mode .. */
      savestatus=set->status;		/* .. and status */
      set->round=rmode;			/* set mode */
      decFloatZero(&zero);		/* make 0E+0 */
      set->status=0;			/* clear */
      decFloatQuantize(&result, df, &zero, set); /* [this may fail] */
      set->round=saveround;		/* restore rounding mode .. */
      if (exact) set->status|=savestatus; /* include Inexact */
       else set->status=savestatus;	/* .. or just original status */
      }
  
    /* only the last four declets of the coefficient can contain */
    /* non-zero; check for others (and also NaN or Infinity from the */
    /* Quantize) first (see DFISZERO for explanation): */
    /* decFloatShow(&result, "sofar"); */
    #if DOUBLE
    if ((DFWORD(&result, 0)&0x1c03ff00)!=0
     || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
    #elif QUAD
    if ((DFWORD(&result, 2)&0xffffff00)!=0
     ||  DFWORD(&result, 1)!=0
     || (DFWORD(&result, 0)&0x1c003fff)!=0
     || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
    #endif
      set->status|=DEC_Invalid_operation; /* Invalid or out of range */
      return 0;
      }
    /* get last twelve digits of the coefficent into hi & ho, base */
    /* 10**9 (see GETCOEFFBILL): */
    sourlo=DFWORD(&result, DECWORDS-1);
    lo=DPD2BIN0[sourlo&0x3ff]
      +DPD2BINK[(sourlo>>10)&0x3ff]
      +DPD2BINM[(sourlo>>20)&0x3ff];
    sourpen=DFWORD(&result, DECWORDS-2);
    hi=DPD2BIN0[((sourpen<<2) | (sourlo>>30))&0x3ff];
  
    /* according to request, check range carefully */
    if (unsign) {
      if (hi>4 || (hi==4 && lo>294967295) || (hi+lo!=0 && DFISSIGNED(&result))) {
        set->status|=DEC_Invalid_operation; /* out of range */
        return 0;
        }
      return hi*BILLION+lo;
      }
    /* signed */
    if (hi>2 || (hi==2 && lo>147483647)) {
      /* handle the usual edge case */
      if (lo==147483648 && hi==2 && DFISSIGNED(&result)) return 0x80000000;
      set->status|=DEC_Invalid_operation; /* truly out of range */
      return 0;
      }
    i=hi*BILLION+lo;
    if (DFISSIGNED(&result)) i=-i;
    return (uInt)i;
    } /* decToInt32 */
  
  /* ------------------------------------------------------------------ */
  /* decToIntegral -- local routine to effect ToIntegral value	      */
  /*								      */
  /*   result gets the result					      */
  /*   df	    is the decFloat to round				      */
  /*   set    is the context					      */
  /*   rmode  is the rounding mode to use				      */
  /*   exact  is 1 if Inexact should be signalled			      */
  /*   returns result						      */
  /* ------------------------------------------------------------------ */
  static decFloat * decToIntegral(decFloat *result, const decFloat *df,
  				decContext *set, enum rounding rmode,
  				Flag exact) {
    Int  exp;			   /* exponent */
    uInt sourhi;			   /* top word from source decFloat */
    enum rounding saveround;	   /* saver */
    uInt savestatus;		   /* .. */
    decFloat zero;		   /* work */
  
    /* Start decoding the argument */
    sourhi=DFWORD(df, 0);		   /* top word */
    exp=DECCOMBEXP[sourhi>>26];	   /* get exponent high bits (in place) */
  
    if (EXPISSPECIAL(exp)) {	   /* is special? */
      /* NaNs are handled as usual */
      if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
      /* must be infinite; return canonical infinity with sign of df */
      return decInfinity(result, df);
      }
  
    /* Here when the argument is finite */
    /* complete extraction of the exponent */
    exp+=GETECON(df)-DECBIAS;		/* .. + continuation and unbias */
  
    if (exp>=0) return decCanonical(result, df); /* already integral */
  
    saveround=set->round;			/* save rounding mode .. */
    savestatus=set->status;		/* .. and status */
    set->round=rmode;			/* set mode */
    decFloatZero(&zero);			/* make 0E+0 */
    decFloatQuantize(result, df, &zero, set); /* 'integrate'; cannot fail */
    set->round=saveround;			/* restore rounding mode .. */
    if (!exact) set->status=savestatus;	/* .. and status, unless exact */
    return result;
    } /* decToIntegral */
  