CbC/CbC_gcc: libdecnumber/decBasic.c comparison

comparison libdecnumber/decBasic.c @ 55:77e2b8dfacca gcc-4.4.5

update it from 4.4.3 to 4.5.0

author	ryoma <e075725@ie.u-ryukyu.ac.jp>
date	Fri, 12 Feb 2010 23:39:51 +0900
parents	a06113de4d67
children	04ced10e8804

comparison

equal deleted inserted replaced

-:c156f1bd5cd9
+:77e2b8dfacca
 /* ------------------------------------------------------------------ */
 /* 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   */
+/* decQuad (but not decSingle).  The main arithmetic operations are   */
-/* here (Add, Subtract, Multiply, FMA, and Division operators).	      */
+/* 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      */
 /* Private constants */
 #define DIVIDE	    0x80000000	   /* Divide operations [as flags] */
 #define REMAINDER   0x40000000	   /* .. */
 #define DIVIDEINT   0x20000000	   /* .. */
-#define REMNEAR	    0x10000000	   /* .. */
+#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 *);
 /* ------------------------------------------------------------------ */
 /* decCanonical -- copy a decFloat, making canonical		      */
 /*								      */
 /*   result gets the canonicalized df				      */
-/*   df	    is the decFloat to copy and make canonical		      */
+/*   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.    */
 /* ------------------------------------------------------------------ */
 inword--;
 encode=DFWORD(result, inword);
 uoff-=32;
 dpd|=encode<<(10-uoff);	   /* get pending bits */
 }
-dpd&=0x3ff;			   /* clear uninteresting 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 */
+encode|=canon<<(uoff-10);    /* insert the canonical form */
 DFWORD(result, inword)=encode;	/* .. and save */
 continue;
 }
 /* straddled words */
 precode=DFWORD(result, inword+1);	/* get previous */
 /* ------------------------------------------------------------------ */
 /* decDivide -- divide operations				      */
 /*								      */
 /*   result gets the result of dividing dfl by dfr:		      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   dfl    is the first decFloat (lhs) 			      */
 /*   dfr    is the second decFloat (rhs)			      */
 /*   set    is the context					      */
-/*   op	    is the operation selector				      */
+/*   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 DIVBASE   ((uInt)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	 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 */
+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;		   /* .. */
+uInt	 *ua, *ud, *uq; 	   /* .. */
-uByte	 *ub;			   /* .. */
+uByte  *ub;			   /* .. */
+uInt	 uiwork;		   /* for macros */
 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
 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 */
+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 */
+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;
 }
 #error Unexpected Quad DIVOPLEN
 #endif
 #endif
 /* set msu and lsu pointers */
-msua=acc+DIVACCLEN-1;	      /* [leading zeros removed below] */
+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 */
+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 */
 	/* 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 */
+	if (*ud==*ua) { 			  /* all compared equal */
 	  *lsuq+=1;				  /* increment result */
 	  msua=lsua;				  /* collapse acc units */
 	  *msua=0;				  /* .. to a zero */
 	  break;
 	  }
 	       + *(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 */
+	    mul/=*msud;       /* no more to the right */
 	    }
 	   else {
-	    mul=(uLong)(*msua) * (uInt)(DIVBASE<<2) + (*(msua-1)<<2);
+	    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);
 /* 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 */
+uInt top, mid, rem; 		/* work */
 if (*uq==0) {			/* no split needed */
-UINTAT(ub)=0;			/* clear 9 BCD8s */
+UBFROMUI(ub, 0);			/* clear 9 BCD8s */
-UINTAT(ub+4)=0;			/* .. */
+UBFROMUI(ub+4, 0);		/* .. */
 *(ub+8)=0;			/* .. */
 continue;
 }
 /* *uq is non-zero -- split the base-billion digit into */
 /* hi, mid, and low three-digits */
 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]);
+UBFROMUI(ub,   UBTOUI(&BIN2BCD8[top*4]));
-UINTAT(ub+3)=UINTAT(&BIN2BCD8[mid*4]);
+UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
-UINTAT(ub+6)=UINTAT(&BIN2BCD8[rem*4]);
+UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
 } /* BCD conversion loop */
-ub--;					/* -> lsu */
+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;
 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 */
+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 (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(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 */
 /* ------------------------------------------------------------------ */
 /* 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)				      */
+/*   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).		      */
 /* 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 */
+/* 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 MULTBASE  ((uInt)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 */
 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;			   /* .. */
+uByte  *ub;			   /* .. */
+uInt	 uiwork;		   /* for macros */
 #if DECUSE64
-uLong	 accl[MULACCLEN];	   /* lazy accumulator (base-billion+) */
+uLong  accl[MULACCLEN];	   /* lazy accumulator (base-billion+) */
-uLong	 *pl;			   /* work -> lazy accumulator */
+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 */
 #if DECUSE64
 /* zero the accumulator */
 #if MULACCLEN==4
 accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0;
-#else					     /* use a loop */
+#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
 /* 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 */
+/*   DOUBLE	 2    29  32  <2*10**18    0.63       1 */
-/*   QUAD	 4    30  31  <4*10**18	   1.17	      2 */
+/*   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. */
 #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) {
+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 */
 /* 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 */
+/* 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) */
 /* 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 */
+/*   DOUBLE	 2    29  32  <2*10**18    0.63       1 */
-/*   QUAD	 4    30  31  <4*10**18	   1.17	      2 */
+/*   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. */
 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 */
+for (pa=acc;; pa++) { 		/* each low uInt */
 uInt hi, lo;			/* words of exact multiply result */
 uInt hop, estlo;			/* work */
 #if QUAD
-uInt esthi;				/* .. */
+uInt esthi; 			/* .. */
 #endif
 lo=*pa;
-hi=*(pa+MULACCLEN);			/* top 32 bits */
+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 */
 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 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]);
+UBFROMUI(ub,   UBTOUI(&BIN2BCD8[top*4]));
-UINTAT(ub+3)=UINTAT(&BIN2BCD8[mid*4]);
+UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
-UINTAT(ub+6)=UINTAT(&BIN2BCD8[rem*4]);
+UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
 }
 else {				/* *pa==0 */
-UINTAT(ub)=0;			/* clear 9 BCD8s */
+UBFROMUI(ub, 0);			/* clear 9 BCD8s */
-UINTAT(ub+4)=0;			/* .. */
+UBFROMUI(ub+4, 0);		/* .. */
 *(ub+8)=0;			/* .. */
 }
 if (pa==acc) break;
 } /* BCD conversion loop */
 /* ------------------------------------------------------------------ */
 /* decFloatAbs -- absolute value, heeding NaNs, etc.		      */
 /*								      */
 /*   result gets the canonicalized df with sign 0		      */
-/*   df	    is the decFloat to abs				      */
+/*   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	      */
 /* ------------------------------------------------------------------ */
 /* decFloatAdd -- add two decFloats				      */
 /*								      */
 /*   result gets the result of adding dfl and dfr:		      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   dfl    is the first decFloat (lhs) 			      */
 /*   dfr    is the second decFloat (rhs)			      */
 /*   set    is the context					      */
 /*   returns result						      */
 /*								      */
 /* ------------------------------------------------------------------ */
+#if QUAD
+/* Table for testing MSDs for fastpath elimination; returns the MSD of */
+/* a decDouble or decQuad (top 6 bits tested) ignoring the sign. */
+/* Infinities return -32 and NaNs return -128 so that summing the two */
+/* MSDs also allows rapid tests for the Specials (see code below). */
+const Int DECTESTMSD[64]={
+0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5,   6,    7,
+0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128,
+0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5,   6,    7,
+0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128};
+#else
+/* The table for testing MSDs is shared between the modules */
+extern const Int DECTESTMSD[64];
+#endif
 decFloat * decFloatAdd(decFloat *result,
 		       const decFloat *dfl, const decFloat *dfr,
 		       decContext *set) {
 bcdnum num;			   /* for final conversion */
-Int	 expl, expr;		   /* left and right exponents */
+Int	 bexpl, bexpr;		   /* left and right biased exponents */
-uInt	 *ui, *uj;		   /* work */
+uByte  *ub, *us, *ut; 	   /* work */
-uByte	 *ub;			   /* .. */
+uInt	 uiwork;		   /* for macros */
+#if QUAD
+uShort uswork;		   /* .. */
+#endif
 uInt sourhil, sourhir;	   /* top words from source decFloats */
-				   /* [valid only until specials */
+				   /* [valid only through end of */
-				   /* handled or exponents decoded] */
+				   /* fastpath code -- before swap] */
 uInt diffsign;		   /* non-zero if signs differ */
 uInt carry;			   /* carry: 0 or 1 before add loop */
-Int  overlap;			   /* coefficient overlap (if full) */
+Int  overlap; 		   /* coefficient overlap (if full) */
+Int  summ;			   /* sum of the MSDs */
 /* 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 */
+#define CARRYPAT 0x01000000    /* carry=1 pattern */
 #else
-#define CARRYPAT 0x00000001	   /* carry=1 pattern */
+#define CARRYPAT 0x00000001    /* carry=1 pattern */
 #endif
 /* Start decoding the arguments */
-/* the initial exponents are placed into the opposite Ints to */
+/* 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) */
+/* whichever path is taken (thereby minimising worst-case path). */
+/* The calculated exponents will be nonsense when the arguments are */
+/* Special, but are not used in that path */
 sourhil=DFWORD(dfl, 0);	   /* LHS top word */
-expr=DECCOMBEXP[sourhil>>26];	   /* get exponent high bits (in place) */
+summ=DECTESTMSD[sourhil>>26];    /* get first MSD for testing */
+bexpr=DECCOMBEXP[sourhil>>26];   /* get exponent high bits (in place) */
+bexpr+=GETECON(dfl);		   /* .. + continuation */
 sourhir=DFWORD(dfr, 0);	   /* RHS top word */
-expl=DECCOMBEXP[sourhir>>26];
+summ+=DECTESTMSD[sourhir>>26];   /* sum MSDs for testing */
+bexpl=DECCOMBEXP[sourhir>>26];
+bexpl+=GETECON(dfr);
+/* here bexpr has biased exponent from lhs, and vice versa */
 diffsign=(sourhil^sourhir)&DECFLOAT_Sign;
-if (EXPISSPECIAL(expl | expr)) { /* either is special? */
+/* now determine whether to take a fast path or the full-function */
-if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
+/* slow path.  The slow path must be taken when: */
-/* one or two infinities */
+/*   -- both numbers are finite, and: */
-/* two infinities with different signs is invalid */
+/*	     the exponents are different, or */
-if (diffsign && DFISINF(dfl) && DFISINF(dfr))
+/*	     the signs are different, or */
-return decInvalid(result, set);
+/*	     the sum of the MSDs is >8 (hence might overflow) */
-if (DFISINF(dfl)) return decInfinity(result, dfl); /* LHS is infinite */
+/* specialness and the sum of the MSDs can be tested at once using */
-return decInfinity(result, dfr);		       /* RHS must be Infinite */
+/* the summ value just calculated, so the test for specials is no */
-}
+/* longer on the worst-case path (as of 3.60) */
-/* Here when both arguments are finite */
+if (summ<=8) {		   /* MSD+MSD is good, or there is a special */
+if (summ<0) {		   /* there is a special */
-/* complete exponent gathering (keeping swapped) */
+/* Inf+Inf would give -64; Inf+finite is -32 or higher */
-expr+=GETECON(dfl)-DECBIAS;	   /* .. + continuation and unbias */
+if (summ<-64) return decNaNs(result, dfl, dfr, set);  /* one or two NaNs */
-expl+=GETECON(dfr)-DECBIAS;
+/* two infinities with different signs is invalid */
-/* here expr has exponent from lhs, and vice versa */
+if (summ==-64 && diffsign) return decInvalid(result, set);
+if (DFISINF(dfl)) return decInfinity(result, dfl);    /* LHS is infinite */
+return decInfinity(result, dfr);			    /* RHS must be Inf */
+}
+/* Here when both arguments are finite; fast path is possible */
+/* (currently only for aligned and same-sign) */
+if (bexpr==bexpl && !diffsign) {
+uInt tac[DECLETS+1];		/* base-1000 coefficient */
+uInt encode;			/* work */
+/* Get one coefficient as base-1000 and add the other */
+GETCOEFFTHOU(dfl, tac);		/* least-significant goes to [0] */
+ADDCOEFFTHOU(dfr, tac);
+/* here the sum of the MSDs (plus any carry) will be <10 due to */
+/* the fastpath test earlier */
+/* construct the result; low word is the same for both formats */
+encode =BIN2DPD[tac[0]];
+encode|=BIN2DPD[tac[1]]<<10;
+encode|=BIN2DPD[tac[2]]<<20;
+encode|=BIN2DPD[tac[3]]<<30;
+DFWORD(result, (DECBYTES/4)-1)=encode;
+/* collect next two declets (all that remains, for Double) */
+encode =BIN2DPD[tac[3]]>>2;
+encode|=BIN2DPD[tac[4]]<<8;
+#if QUAD
+/* complete and lay out middling words */
+encode|=BIN2DPD[tac[5]]<<18;
+encode|=BIN2DPD[tac[6]]<<28;
+DFWORD(result, 2)=encode;
+encode =BIN2DPD[tac[6]]>>4;
+encode|=BIN2DPD[tac[7]]<<6;
+encode|=BIN2DPD[tac[8]]<<16;
+encode|=BIN2DPD[tac[9]]<<26;
+DFWORD(result, 1)=encode;
+/* and final two declets */
+encode =BIN2DPD[tac[9]]>>6;
+encode|=BIN2DPD[tac[10]]<<4;
+#endif
+/* add exponent continuation and sign (from either argument) */
+encode|=sourhil & (ECONMASK | DECFLOAT_Sign);
+/* create lookup index = MSD + top two bits of biased exponent <<4 */
+tac[DECLETS]|=(bexpl>>DECECONL)<<4;
+encode|=DECCOMBFROM[tac[DECLETS]]; /* add constructed combination field */
+DFWORD(result, 0)=encode; 	 /* complete */
+/* decFloatShow(result, ">"); */
+return result;
+} /* fast path OK */
+/* drop through to slow path */
+} /* low sum or Special(s) */
+/* Slow path required -- arguments are finite and might overflow,   */
+/* or require alignment, or might have different signs	      */
 /* now swap either exponents or argument pointers */
-if (expl<=expr) {
+if (bexpl<=bexpr) {
 /* original left is bigger */
-Int expswap=expl;
+Int bexpswap=bexpl;
-expl=expr;
+bexpl=bexpr;
-expr=expswap;
+bexpr=bexpswap;
 /* 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, */
+/* [here dfl and bexpl 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 */
 /* 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 */
+UBFROMUS(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.exponent=bexpl-DECBIAS;
 tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
 decShowNum(&tum, "dflx");}
 #endif
 /* if signs differ, take ten's complement of lhs (here the */
 /* 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);
+UBFROMUI(acc+ 4, 0x09090909-UBTOUI(acc+ 4));
-UINTAT(acc+ 8)=0x09090909-UINTAT(acc+ 8);
+UBFROMUI(acc+ 8, 0x09090909-UBTOUI(acc+ 8));
-UINTAT(acc+12)=0x09090909-UINTAT(acc+12);
+UBFROMUI(acc+12, 0x09090909-UBTOUI(acc+12));
-UINTAT(acc+16)=0x09090909-UINTAT(acc+16);
+UBFROMUI(acc+16, 0x09090909-UBTOUI(acc+16));
 #if QUAD
-UINTAT(acc+20)=0x09090909-UINTAT(acc+20);
+UBFROMUI(acc+20, 0x09090909-UBTOUI(acc+20));
-UINTAT(acc+24)=0x09090909-UINTAT(acc+24);
+UBFROMUI(acc+24, 0x09090909-UBTOUI(acc+24));
-UINTAT(acc+28)=0x09090909-UINTAT(acc+28);
+UBFROMUI(acc+28, 0x09090909-UBTOUI(acc+28));
-UINTAT(acc+32)=0x09090909-UINTAT(acc+32);
+UBFROMUI(acc+32, 0x09090909-UBTOUI(acc+32));
-UINTAT(acc+36)=0x09090909-UINTAT(acc+36);
+UBFROMUI(acc+36, 0x09090909-UBTOUI(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);
+overlap=DECPMAX-(bexpl-bexpr);
 #if DECTRACE
-printf("exps: %ld %ld\n", (LI)expl, (LI)expr);
+printf("exps: %ld %ld\n", (LI)(bexpl-DECBIAS), (LI)(bexpr-DECBIAS));
 printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry);
 #endif
 if (overlap<=0) {			/* no overlap possible */
 uInt gap;				/* local work */
 /* 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;
+bexpr+=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 */
+ut=acc+COFF+DECPMAX;		/* start of gap */
-for (; ui<&UINTAT(ub); ui++) *ui=0; /* mind the gap */
+for (; ut<ub; ut+=4) UBFROMUI(ut, 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 */
 } /* 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 */
+if (overlap==DECPMAX) {		/* aligned */
-/* Fill the prefix gap with 0s; 8 will cover most common */
+ub=buf+COFF;			/* where msd will go */
-/* unalignments, so start with direct assignments (a loop is */
+#if QUAD
-/* then used for any remaining -- the loop (and the one in a */
+UBFROMUS(buf+4, 0);		/* clear quad's 00 */
-/* moment) is not then on the critical path because the number */
+#endif
-/* of additions is reduced by (at least) two in this case) */
+GETCOEFF(dfr, ub);		/* decode from decFloat */
-UINTAT(buf+4)=0;			/* [clears decQuad 00 too] */
+}
-UINTAT(buf+8)=0;
+else {				/* unaligned */
-if (ub>buf+12) {
+ub=buf+COFF+DECPMAX-overlap;	/* where MSD will go */
-ui=&UINTAT(buf+12);		/* start of any remaining */
+/* Fill the prefix gap with 0s; 8 will cover most common */
-for (; ui<&UINTAT(ub); ui++) *ui=0; /* fill them */
+/* unalignments, so start with direct assignments (a loop is */
-}
+/* then used for any remaining -- the loop (and the one in a */
-GETCOEFF(dfr, ub);			/* decode from decFloat */
+/* moment) is not then on the critical path because the number */
+/* of additions is reduced by (at least) two in this case) */
-/* now move tail of rhs across to main acc; again use direct */
+UBFROMUI(buf+4, 0);		/* [clears decQuad 00 too] */
-/* assignment for 8 digits-worth */
+UBFROMUI(buf+8, 0);
-UINTAT(acc+COFF+DECPMAX)=UINTAT(buf+COFF+DECPMAX);
+if (ub>buf+12) {
-UINTAT(acc+COFF+DECPMAX+4)=UINTAT(buf+COFF+DECPMAX+4);
+	ut=buf+12;			/* start any remaining */
-if (buf+COFF+DECPMAX+8<ub+DECPMAX) {
+	for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* fill them */
-uj=&UINTAT(buf+COFF+DECPMAX+8);	/* source */
+	}
-ui=&UINTAT(acc+COFF+DECPMAX+8);	/* target */
+GETCOEFF(dfr, ub);		/* decode from decFloat */
-for (; uj<&UINTAT(ub+DECPMAX); ui++, uj++) *ui=*uj;
-}
+/* now move tail of rhs across to main acc; again use direct */
+/* copies for 8 digits-worth */
+UBFROMUI(acc+COFF+DECPMAX,   UBTOUI(buf+COFF+DECPMAX));
+UBFROMUI(acc+COFF+DECPMAX+4, UBTOUI(buf+COFF+DECPMAX+4));
+if (buf+COFF+DECPMAX+8<ub+DECPMAX) {
+	us=buf+COFF+DECPMAX+8;		/* source */
+	ut=acc+COFF+DECPMAX+8;		/* target */
+	for (; us<ub+DECPMAX; us+=4, ut+=4) UBFROMUI(ut, UBTOUI(us));
+	}
+} /* unaligned */
 ulsd=acc+(ub-buf+DECPMAX-1);	/* update LSD pointer */
-/* now do the add of the non-tail; this is all nicely aligned, */
+/* 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 */
+/* zero 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] */
+/* [byte-by-byte add, here, is about 15% slower total effect 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) */
+ut=acc+COFF+DECPMAX-4;		/* target LSW (acc) */
-uj=&UINTAT(buf+COFF+DECPMAX-4);	/* source LSW (buf, to add to acc) */
+us=buf+COFF+DECPMAX-4;		/* source LSW (buf, to add to acc) */
 #if !DECLITEND
-for (; ui>=&UINTAT(acc+4); ui--, uj--) {
+for (; ut>=acc+4; ut-=4, us-=4) {	/* big-endian add loop */
 /* bcd8 add */
-carry+=*uj;			/* rhs + carry */
+carry+=UBTOUI(us);		/* rhs + carry */
 if (carry==0) continue;		/* no-op */
-carry+=*ui;			/* lhs */
+carry+=UBTOUI(ut);		/* lhs */
 /* Big-endian BCD adjust (uses internal carry) */
 carry+=0x76f6f6f6;		/* note top nibble not all bits */
-*ui=(carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4); /* BCD adjust */
+/* apply BCD adjust and save */
+UBFROMUI(ut, (carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4));
 carry>>=31;			/* true carry was at far left */
 } /* add loop */
 #else
-for (; ui>=&UINTAT(acc+4); ui--, uj--) {
+for (; ut>=acc+4; ut-=4, us-=4) {	/* little-endian add loop */
 /* bcd8 add */
-carry+=*uj;			/* rhs + carry */
+carry+=UBTOUI(us);		/* rhs + carry */
 if (carry==0) continue;		/* no-op [common if unaligned] */
-carry+=*ui;			/* lhs */
+carry+=UBTOUI(ut);		/* 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 */
+UBFROMUI(ut, 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;
 /* if big-endian */
 #if !DECLITEND
 *(ulsd+1)=0;
 #endif
 /* there are always at least four coefficient words */
-UINTAT(umsd)   =0x09090909-UINTAT(umsd);
+UBFROMUI(umsd,	0x09090909-UBTOUI(umsd));
-UINTAT(umsd+4) =0x09090909-UINTAT(umsd+4);
+UBFROMUI(umsd+4,	0x09090909-UBTOUI(umsd+4));
-UINTAT(umsd+8) =0x09090909-UINTAT(umsd+8);
+UBFROMUI(umsd+8,	0x09090909-UBTOUI(umsd+8));
-UINTAT(umsd+12)=0x09090909-UINTAT(umsd+12);
+UBFROMUI(umsd+12, 0x09090909-UBTOUI(umsd+12));
 #if DOUBLE
 	#define BNEXT 16
 #elif QUAD
-	UINTAT(umsd+16)=0x09090909-UINTAT(umsd+16);
+	UBFROMUI(umsd+16, 0x09090909-UBTOUI(umsd+16));
-	UINTAT(umsd+20)=0x09090909-UINTAT(umsd+20);
+	UBFROMUI(umsd+20, 0x09090909-UBTOUI(umsd+20));
-	UINTAT(umsd+24)=0x09090909-UINTAT(umsd+24);
+	UBFROMUI(umsd+24, 0x09090909-UBTOUI(umsd+24));
-	UINTAT(umsd+28)=0x09090909-UINTAT(umsd+28);
+	UBFROMUI(umsd+28, 0x09090909-UBTOUI(umsd+28));
-	UINTAT(umsd+32)=0x09090909-UINTAT(umsd+32);
+	UBFROMUI(umsd+32, 0x09090909-UBTOUI(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);
+	UBFROMUI(umsd+BNEXT,   0x09090909-UBTOUI(umsd+BNEXT));
-	UINTAT(umsd+BNEXT+4)=0x09090909-UINTAT(umsd+BNEXT+4);
+	UBFROMUI(umsd+BNEXT+4, 0x09090909-UBTOUI(umsd+BNEXT+4));
 	#if DOUBLE
 	#define BNEXTY (BNEXT+8)
 	#elif QUAD
-	UINTAT(umsd+BNEXT+8)=0x09090909-UINTAT(umsd+BNEXT+8);
+	UBFROMUI(umsd+BNEXT+8,	0x09090909-UBTOUI(umsd+BNEXT+8));
-	UINTAT(umsd+BNEXT+12)=0x09090909-UINTAT(umsd+BNEXT+12);
+	UBFROMUI(umsd+BNEXT+12, 0x09090909-UBTOUI(umsd+BNEXT+12));
 	#define BNEXTY (BNEXT+16)
 	#endif
 	if (ulsd>=umsd+BNEXTY) {	/* very unaligned */
-	  ui=&UINTAT(umsd+BNEXTY);	/* -> continue */
+	  ut=umsd+BNEXTY;		/* -> continue */
-	  for (;;ui++) {
+	  for (;;ut+=4) {
-	    *ui=0x09090909-*ui;		/* invert four digits */
+	    UBFROMUI(ut, 0x09090909-UBTOUI(ut)); /* invert four digits */
-	    if (ui>=&UINTAT(ulsd-3)) break; /* all done */
+	    if (ut>=ulsd-3) break;	/* all done */
 	    }
 	  }
 	}
 /* complete the ten's complement by adding 1 */
 for (ub=ulsd; *ub==9; ub--) *ub=0;
 /* 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 (; UBTOUI(umsd)==0 && umsd+3<ulsd;) umsd+=4;
 	  for (; *umsd==0 && umsd<ulsd;) umsd++;
 	  }
-	if (*umsd==0) {		   /* must be true zero (and diffsign) */
+	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] */
 umsd=acc+3;
 }
 #endif
 } /* same sign */
-num.msd=umsd;			   /* set MSD .. */
+num.msd=umsd; 		   /* set MSD .. */
-num.lsd=ulsd;			   /* .. and LSD */
+num.lsd=ulsd; 		   /* .. and LSD */
-num.exponent=expr;		   /* set exponent to smaller */
+num.exponent=bexpr-DECBIAS;	   /* set exponent to smaller, unbiassed */
 #if DECTRACE
 decFloatShow(dfl, "dfl");
 decFloatShow(dfr, "dfr");
 decShowNum(&num, "postadd");
 /* ------------------------------------------------------------------ */
 /* decFloatAnd -- logical digitwise AND of two decFloats	      */
 /*								      */
 /*   result gets the result of ANDing dfl and dfr		      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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	      */
+/* 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) {
 /* ------------------------------------------------------------------ */
 /* decFloatCanonical -- copy a decFloat, making canonical	      */
 /*								      */
 /*   result gets the canonicalized df				      */
-/*   df	    is the decFloat to copy and make canonical		      */
+/*   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) {
 } /* decFloatCanonical */
 /* ------------------------------------------------------------------ */
 /* decFloatClass -- return the class of a decFloat		      */
 /*								      */
-/*   df is the decFloat to test					      */
+/*   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)) {
 } /* decFloatClass */
 /* ------------------------------------------------------------------ */
 /* decFloatClassString -- return the class of a decFloat as a string  */
 /*								      */
-/*   df is the decFloat to test					      */
+/*   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_SNAN)	       return DEC_ClassString_SN;
 return DEC_ClassString_UN;	       /* Unknown */
 } /* decFloatClassString */
 /* ------------------------------------------------------------------ */
-/* decFloatCompare -- compare two decFloats; quiet NaNs allowed	      */
+/* decFloatCompare -- compare two decFloats; quiet NaNs allowed       */
 /*								      */
 /*   result gets the result of comparing dfl and dfr		      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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,
 /* ------------------------------------------------------------------ */
 /* decFloatCompareSignal -- compare two decFloats; all NaNs signal    */
 /*								      */
 /*   result gets the result of comparing dfl and dfr		      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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,
 /* ------------------------------------------------------------------ */
 /* decFloatCompareTotal -- compare two decFloats with total ordering  */
 /*								      */
 /*   result gets the result of comparing dfl and dfr		      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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 */
+Int  comp;				     /* work */
+uInt uiwork;				     /* for macros */
+#if QUAD
+uShort uswork;			     /* .. */
+#endif
 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 */
+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 */
+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 */
+	UBFROMUS(bufl, 0);
-	USHORTAT(ub)=0;
+	UBFROMUS(bufr, 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 */
+	uInt ui=UBTOUI(ub);
+	if (ui==UBTOUI(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;		     /* .. */
 /* ------------------------------------------------------------------ */
 /* decFloatCompareTotalMag -- compare magnitudes with total ordering  */
 /*								      */
 /*   result gets the result of comparing abs(dfl) and abs(dfr)	      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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) {
 } /* decFloatCopyAbs */
 /* ------------------------------------------------------------------ */
 /* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
 /*								      */
-/*   result gets the copy of dfl with sign bit inverted		      */
+/*   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. */
 /* ------------------------------------------------------------------ */
 DFBYTE(result, 0)^=0x80;		/* invert sign bit */
 return result;
 } /* decFloatCopyNegate */
 /* ------------------------------------------------------------------ */
-/* decFloatCopySign -- copy a decFloat with the sign of another	      */
+/* decFloatCopySign -- copy a decFloat with the sign of another       */
 /*								      */
-/*   result gets the result of copying dfl with the sign of dfr	      */
+/*   result gets the result of copying dfl with the sign of dfr       */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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. */
 /* ------------------------------------------------------------------ */
 #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;				  \
+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 */
 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 */
+dpdlenchk(1, sourlo>>30);  /* very rare */
 dpdlendun(2, sourlo>>20);
 } /* [cannot drop through] */
 dpdlenchk(3, sourlo>>10);
 dpdlendun(4, sourlo);
 /* [cannot drop through] */
 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 */
+dpdlenchk(7, sourlo>>30);  /* very rare */
 dpdlendun(8, sourlo>>20);
 } /* [cannot drop through] */
 dpdlenchk(9, sourlo>>10);
 dpdlendun(10, sourlo);
 /* [cannot drop through] */
 /* ------------------------------------------------------------------ */
 /* decFloatDivide -- divide a decFloat by another		      */
 /*								      */
 /*   result gets the result of dividing dfl by dfr:		      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   dfl    is the first decFloat (lhs) 			      */
 /*   dfr    is the second decFloat (rhs)			      */
 /*   set    is the context					      */
 /*   returns result						      */
 /*								      */
 /* ------------------------------------------------------------------ */
 /* ------------------------------------------------------------------ */
 /* decFloatDivideInteger -- integer divide a decFloat by another      */
 /*								      */
 /*   result gets the result of dividing dfl by dfr:		      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   dfl    is the first decFloat (lhs) 			      */
 /*   dfr    is the second decFloat (rhs)			      */
 /*   set    is the context					      */
 /*   returns result						      */
 /*								      */
 /* ------------------------------------------------------------------ */
 /* ------------------------------------------------------------------ */
 /* 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)				      */
+/*   dfl    is the first decFloat (lhs) 			      */
 /*   dfr    is the second decFloat (rhs)			      */
-/*   dff    is the final decFloat (fhs)				      */
+/*   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)
+#define FMALEN (ROUNDUP4(1+ (DECPMAX9*18+1) +DECPMAX+2))
-uByte	 acc[FMALEN];		   /* for multiplied coefficient in BCD */
+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 */
+uByte  coe[ROUNDUP4(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 */
+uInt	 hipad; 		   /* pad digit for hi if needed */
 Int	 padding;		   /* excess exponent */
-uInt	 carry;			   /* +1 for ten's complement and during add */
+uInt	 carry; 		   /* +1 for ten's complement and during add */
-uByte	 *ub, *uh, *ul;		   /* work */
+uByte  *ub, *uh, *ul; 	   /* work */
+uInt	 uiwork;		   /* for macros */
 /* 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 */
 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 */
+/* has the higher exponent and lo points to the other [don't care, */
-/* the same] */
+/* if the same].  One coefficient will be in acc, the other in coe. */
 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 */
+/* and make testing for zero trivial (tests are safe because acc */
-for (; UINTAT(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4;
+/* and coe are rounded up to uInts) */
+for (; UBTOUI(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 (; UBTOUI(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 */
+if (*hi->msd==0) {			     /* 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 */
+if (*lo->msd==0) {		     /* 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 lo is zero then padding with zeros may be needed, below) */
 /* if signs differ, take the ten's complement of hi (zeros to the */
-/* right do not matter because the complement of zero is zero); */
+/* right do not matter because the complement of zero is zero); the */
-/* the +1 is done later, as part of the addition, inserted at 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->msd; uh<hi->lsd-3; uh+=4) UBFROMUI(uh, 0x09090909-UBTOUI(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 */
 /* 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 */
+/* starting from the far right; this could be either hi or lo, and */
+/* will be aligned */
 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 hilen=(Int)(hi->lsd-hi->msd+1); /* length of hi */
-Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
+Int lolen=(Int)(lo->lsd-lo->msd+1); /* and of lo */
-Int reduce=newexp-lo->exponent;
-if (reduce>0) {			/* [= case gives reduce=0 nop] */
+if (hilen+padding-lolen > DECPMAX+2) {   /* can reduce lo to single */
+/* make sure it is virtually at least DECPMAX from hi->msd, at */
+/* least to right of hi->lsd (in case of destructive subtract), */
+/* and separated by at least two digits from either of those */
+/* (the 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 to: */
+/*   hi:	 1000000000000000000		     E-16 */
+/*   lo:	  .................1000000000000000  E-16 */
+Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
 /* printf("FMA reduce: %ld\n", (LI)reduce); */
-if (reduce>=lolen) {		/* eating all */
+lo->lsd=lo->msd;			     /* to single digit [maybe 0] */
-	lo->lsd=lo->msd;		/* reduce to single digit */
+lo->exponent=newexp;		     /* new lowest exponent */
-	lo->exponent=newexp;		/* [known to be non-zero] */
+padding=hi->exponent-lo->exponent;     /* recalculate */
-	}
+ul=lo->lsd;			     /* .. and repoint */
-else { /* < */
+}
-	uByte *up=lo->lsd;
-	lo->lsd=lo->lsd-reduce;
+/* padding is still > 0, but will fit in acc (less leading carry slot) */
-	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);
+if (hilen+padding+1>FMALEN)
+	printf("FMA excess hilen+padding: %ld+%ld \n", (LI)hilen, (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-3>=lo->msd && padding>3; padding-=4, ul-=4, ub-=4) {
+UBFROMUI(ub-3, UBTOUI(ul-3));	     /* [cannot overlap] */
+}
 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... */
+/* dow do the add from hi->lsd to the left */
+/* [bytewise, because either operand can run out at any time] */
+/* carry was set up depending on ten's complement above */
+/* first assume both operands have some digits */
 for (;; ub--) {
-uInt hid, lod;
+if (uh<hi->msd || ul<lo->msd) break;
-if (uh<hi->msd) {
+*ub=(uByte)(carry+(*uh--)+(*ul--));
-if (ul<lo->msd) break;
+carry=0;
-hid=hipad;
+if (*ub<10) continue;
-}
+*ub-=10;
-else hid=*uh--;
+carry=1;
-if (ul<lo->msd) lod=0;
+} /* both loop */
-else lod=*ul--;
-*ub=(uByte)(carry+hid+lod);
+if (ul<lo->msd) {	      /* to left of lo */
-if (*ub<10) carry=0;
+for (;; ub--) {
-else {
+if (uh<hi->msd) break;
+*ub=(uByte)(carry+(*uh--));  /* [+0] */
+carry=0;
+if (*ub<10) continue;
 *ub-=10;
 carry=1;
-}
+} /* hi loop */
-} /* addition loop */
+}
+else {		      /* to left of hi */
+for (;; ub--) {
+if (ul<lo->msd) break;
+*ub=(uByte)(carry+hipad+(*ul--));
+carry=0;
+if (*ub<10) continue;
+*ub-=10;
+carry=1;
+} /* lo 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;
 } /* 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->msd; ul<lo->lsd-3; ul+=4) UBFROMUI(ul, 0x09090909-UBTOUI(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 (; UBTOUI(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 */
+#if DECCHECK
+/* assert no left underrun */
+if (lo->msd<acc) {
+printf("FMA underrun by %ld \n", (LI)(acc-lo->msd));
+}
+#endif
 return decFinalize(result, lo, set);	/* round, check, and lay out */
 } /* decFloatFMA */
 /* ------------------------------------------------------------------ */
-/* decFloatFromInt -- initialise a decFloat from an Int		      */
+/* decFloatFromInt -- initialise a decFloat from an Int 	      */
 /*								      */
 /*   result gets the converted Int				      */
 /*   n	    is the Int to convert				      */
 /*   returns result						      */
 /*								      */
 /* ------------------------------------------------------------------ */
 /* decFloatInvert -- logical digitwise INVERT of a decFloat	      */
 /*								      */
 /*   result gets the result of INVERTing df			      */
-/*   df	    is the decFloat to invert				      */
+/*   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.   */
 } /* decFloatInvert */
 /* ------------------------------------------------------------------ */
 /* decFloatIs -- decFloat tests (IsSigned, etc.)		      */
 /*								      */
-/*   df is the decFloat to test					      */
+/*   df is the decFloat to test 				      */
-/*   returns 0 or 1 in an int32_t				      */
+/*   returns 0 or 1 in a uInt					      */
 /*								      */
 /* 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  */
+/* is a little cleaner (and they can be referred to here by the       */
-/* names)							      */
+/* generic names)						      */
 /* ------------------------------------------------------------------ */
 uInt decFloatIsCanonical(const decFloat *df) {
 if (DFISSPECIAL(df)) {
 if (DFISINF(df)) {
 if (DFWORD(df, 0)&ECONMASK) return 0;  /* exponent continuation */
 uInt decFloatIsZero(const decFloat *df) {
 return DFISZERO(df);
 } /* decFloatIs... */
 /* ------------------------------------------------------------------ */
-/* decFloatLogB -- return adjusted exponent, by 754r rules	      */
+/* decFloatLogB -- return adjusted exponent, by 754 rules	      */
 /*								      */
 /*   result gets the adjusted exponent as an integer, or a NaN etc.   */
-/*   df	    is the decFloat to be examined			      */
+/*   df     is the decFloat to be examined			      */
 /*   set    is the context					      */
 /*   returns result						      */
 /*								      */
 /* Notable cases:						      */
 /*   A<0 -> Use |A|						      */
 			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 */
+return decInfinity(result, result);      /* canonical +Infinity */
 }
 if (DFISZERO(df)) {
-set->status|=DEC_Division_by_zero;	     /* as per 754r */
+set->status|=DEC_Division_by_zero;	     /* as per 754 */
 DFWORD(result, 0)=DECFLOAT_Sign;	     /* make negative */
-return decInfinity(result, result);	     /* canonical -Infinity */
+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 */
 #endif
 return result;
 } /* decFloatLogB */
 /* ------------------------------------------------------------------ */
-/* decFloatMax -- return maxnum of two operands			      */
+/* decFloatMax -- return maxnum of two operands 		      */
 /*								      */
 /*   result gets the chosen decFloat				      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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.		      */
 /* ------------------------------------------------------------------ */
 /* decFloatMaxMag -- return maxnummag of two operands		      */
 /*								      */
 /*   result gets the chosen decFloat				      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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 */
 if (comp<0) return decCanonical(result, dfr);
 return decFloatMax(result, dfl, dfr, set);
 } /* decFloatMaxMag */
 /* ------------------------------------------------------------------ */
-/* decFloatMin -- return minnum of two operands			      */
+/* decFloatMin -- return minnum of two operands 		      */
 /*								      */
 /*   result gets the chosen decFloat				      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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.		      */
 /* ------------------------------------------------------------------ */
 /* decFloatMinMag -- return minnummag of two operands		      */
 /*								      */
 /*   result gets the chosen decFloat				      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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 */
 } /* decFloatMinMag */
 /* ------------------------------------------------------------------ */
 /* decFloatMinus -- negate value, heeding NaNs, etc.		      */
 /*								      */
-/*   result gets the canonicalized 0-df				      */
+/*   result gets the canonicalized 0-df 			      */
-/*   df	    is the decFloat to minus				      */
+/*   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).			      */
 } /* decFloatMinus */
 /* ------------------------------------------------------------------ */
 /* decFloatMultiply -- multiply two decFloats			      */
 /*								      */
-/*   result gets the result of multiplying dfl and dfr:		      */
+/*   result gets the result of multiplying dfl and dfr: 	      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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 */
+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 */
 /*   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   */
+/* This is 754 nextdown; Invalid is the only status possible (from    */
 /* an sNaN).							      */
 /* ------------------------------------------------------------------ */
 decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
 			     decContext *set) {
 decFloat delta;			/* tiny increment */
 /* 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 */
+decFloatZero(&delta); 		/* set up tiny delta */
-DFWORD(&delta, DECWORDS-1)=1;		/* coefficient=1 */
+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 */
+saveround=set->round; 		/* save mode */
 set->round=DEC_ROUND_FLOOR;		/* .. round towards -Infinity */
-savestat=set->status;			/* save status */
+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 */
+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     */
+/* This is 754 nextup; Invalid is the only status possible (from      */
 /* an sNaN).							      */
 /* ------------------------------------------------------------------ */
 decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
 			    decContext *set) {
 uInt savestat;			/* saves status */
 /* 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 */
+decFloatZero(&delta); 		/* set up tiny delta */
-DFWORD(&delta, DECWORDS-1)=1;		/* coefficient=1 */
+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 */
+saveround=set->round; 		/* save mode */
-set->round=DEC_ROUND_CEILING;		/* .. round towards +Infinity */
+set->round=DEC_ROUND_CEILING; 	/* .. round towards +Infinity */
-savestat=set->status;			/* save status */
+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 */
+set->round=saveround; 		/* .. and mode */
 return result;
 } /* decFloatNextPlus */
 /* ------------------------------------------------------------------ */
 /* decFloatNextToward -- next towards a 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   */
+/* This is 754-1985 nextafter, as modified during revision (dropped   */
-/* normal number.						      */
+/* from 754-2008); 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 */
 DFWORD(result, 0)|=DECFLOAT_Sign;
 return result;
 }
 saveround=set->round;		     /* save mode */
 set->round=DEC_ROUND_CEILING;	     /* .. round towards +Infinity */
-deltatop=0;				     /* positive delta */
+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 */
+set->round=DEC_ROUND_FLOOR; 	     /* .. round towards -Infinity */
 deltatop=DECFLOAT_Sign;		     /* negative delta */
 }
-savestat=set->status;			     /* save status */
+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 */
+decFloatZero(&delta); 		/* set up tiny delta */
-DFWORD(&delta, DECWORDS-1)=1;		/* coefficient=1 */
+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 */
+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)				      */
+/*   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	      */
+/* 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) {
 } /* decFloatOr */
 /* ------------------------------------------------------------------ */
 /* decFloatPlus -- add value to 0, heeding NaNs, etc.		      */
 /*								      */
-/*   result gets the canonicalized 0+df				      */
+/*   result gets the canonicalized 0+df 			      */
-/*   df	    is the decFloat to plus				      */
+/*   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	      */
+/* 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) {
 /* ------------------------------------------------------------------ */
 /* decFloatQuantize -- quantize a decFloat			      */
 /*								      */
 /*   result gets the result of quantizing dfl to match dfr	      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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   */
 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, *uc;	      /* work */
-uByte *ub;		      /* .. */
 Int	drop;		      /* .. */
 uInt	dpd;		      /* .. */
-uInt	encode;		      /* encoding accumulator */
+uInt	encode; 	      /* encoding accumulator */
 uInt	sourhil, sourhir;     /* top words from source decFloats */
+uInt	uiwork; 	      /* for macros */
+#if QUAD
+uShort uswork;	      /* .. */
+#endif
 /* the following buffer holds the coefficient for manipulation */
-uByte buf[4+DECPMAX*3];     /* + space for zeros to left or right */
+uByte buf[4+DECPMAX*3+2*QUAD];   /* + space for zeros to left or right */
 #if DECTRACE
-bcdnum num;		      /* for trace displays */
+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) */
 num.exponent=explb-DECBIAS;
 num.sign=sourhil & DECFLOAT_Sign;
 decShowNum(&num, "dfl");
 #endif
-if (drop>0) {				/* [most common case] */
+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;
+UBFROMUI(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 */
+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 */
 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;
+	for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(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] */
 /* 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 */
+UBFROMUI(BUFOFF-8, 0);		     /* must be at least 5 */
-for (ui=&UINTAT(BUFOFF-12); ui>&UINTAT(ulsd-DECPMAX-3); ui--) *ui=0;
+for (uc=BUFOFF-12; uc>ulsd-DECPMAX-3; uc-=4) UBFROMUI(uc, 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 */
 #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++) {
+/* note that here zeros to the right are added by fours, so in */
-	*ui=0;
+/* the Quad case this could write 36 zeros if the coefficient has */
-	if (UINTAT(&UBYTEAT(ui)-DECPMAX)!=0) { /* could be bad */
+/* fewer than three significant digits (hence the +2*QUAD for buf) */
+for (uc=BUFOFF+DECPMAX;; uc+=4) {
+	UBFROMUI(uc, 0);
+	if (UBTOUI(uc-DECPMAX)!=0) {		  /* could be bad */
 	  /* if all four digits should be zero, definitely bad */
-	  if (ui<=&UINTAT(BUFOFF+DECPMAX+(-drop)-4))
+	  if (uc<=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)
+	  if ((UBTOUI(uc-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 */
+	if (uc>=BUFOFF+DECPMAX+(-drop)-4) break;  /* done */
 	}
 ulsd=BUFOFF+DECPMAX+(-drop)-1;
 } /* pad and check leading zeros */
 } /* drop<0 */
 /* ------------------------------------------------------------------ */
 /* decFloatReduce -- reduce finite coefficient to minimum length      */
 /*								      */
 /*   result gets the reduced decFloat				      */
-/*   df	    is the source 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.	      */
 /* ------------------------------------------------------------------ */
 /* decFloatRemainder -- integer divide and return remainder	      */
 /*								      */
 /*   result gets the remainder of dividing dfl by dfr:		      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   dfl    is the first decFloat (lhs) 			      */
 /*   dfr    is the second decFloat (rhs)			      */
 /*   set    is the context					      */
 /*   returns result						      */
 /*								      */
 /* ------------------------------------------------------------------ */
 /* ------------------------------------------------------------------ */
 /* decFloatRemainderNear -- integer divide to nearest and remainder   */
 /*								      */
 /*   result gets the remainder of dividing dfl by dfr:		      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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.     */
 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 */
+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 */
 }
 /* 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 */
+savestat=set->status; 		/* record */
 decFinalize(result, &num, set);
-set->status=savestat;			/* restore */
+set->status=savestat; 		/* restore */
 return result;
 } /* decFloatRotate */
 /* ------------------------------------------------------------------ */
 /* decFloatSameQuantum -- test decFloats for same quantum	      */
 /*								      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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.			      */
 /* ------------------------------------------------------------------ */
 if (GETEXP(dfl)==GETEXP(dfr)) return 1; /* biased exponents match */
 return 0;
 } /* decFloatSameQuantum */
 /* ------------------------------------------------------------------ */
-/* decFloatScaleB -- multiply by a power of 10, as per 754r	      */
+/* decFloatScaleB -- multiply by a power of 10, as per 754	      */
 /*								      */
 /*   result gets the result of the operation			      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   dfl    is the first decFloat (lhs) 			      */
-/*   dfr    is the second decFloat (rhs), am integer (with q=0)	      */
+/*   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.	      */
 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 */
+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 */
+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 */
+*result=*dfl; 			     /* copy to target */
 return decFloatSetExponent(result, set, GETEXPUN(result)+expr);
 } /* decFloatScaleB */
 /* ------------------------------------------------------------------ */
 /* decFloatShift -- shift the coefficient of a decFloat left or right */
 /* 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 */
+Int	 shift; 			/* dfr as an Int */
-uByte buf[DECPMAX*2];			/* coefficient + padding */
+uByte  buf[DECPMAX*2];		/* coefficient + padding */
-uInt digits, savestat;		/* work */
+uInt	 digits, savestat;		/* work */
 bcdnum num;				/* .. */
+uInt	 uiwork;			/* for macros */
 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 */
+if (digits>2) return decInvalid(result, set);   /* definitely out of range */
-shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff];	  /* is in bottom declet */
+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 */
+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;
 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 */
+UBFROMUI(buf+DECPMAX, 0);		/* 8 will handle most cases */
-UINTAT(buf+DECPMAX+4)=0;		/* .. */
+UBFROMUI(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 */
+savestat=set->status; 		/* record */
 decFinalize(result, &num, set);
-set->status=savestat;			/* restore */
+set->status=savestat; 		/* restore */
 return result;
 } /* decFloatShift */
 /* ------------------------------------------------------------------ */
-/* decFloatSubtract -- subtract a decFloat from another		      */
+/* decFloatSubtract -- subtract a decFloat from another 	      */
 /*								      */
 /*   result gets the result of subtracting dfr from dfl:	      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   dfl    is the first decFloat (lhs) 			      */
 /*   dfr    is the second decFloat (rhs)			      */
 /*   set    is the context					      */
 /*   returns result						      */
 /*								      */
 /* ------------------------------------------------------------------ */
 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)	      */
+/* decFloatToInt -- round to 32-bit binary integer (4 flavours)       */
 /*								      */
-/*   df	   is the decFloat to round				      */
+/*   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 */
 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)	      */
+/* decFloatToIntegral -- round to integral value (two flavours)       */
 /*								      */
 /*   result gets the result					      */
-/*   df	    is the decFloat to round				      */
+/*   df     is the decFloat to round				      */
 /*   set    is the context					      */
-/*   round  is the rounding mode to use				      */
+/*   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.				      */
 /* ------------------------------------------------------------------ */
 /* ------------------------------------------------------------------ */
 /* decFloatXor -- logical digitwise XOR of two decFloats	      */
 /*								      */
 /*   result gets the result of XORing dfl and dfr		      */
-/*   dfl    is the first decFloat (lhs)				      */
+/*   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	      */
+/* 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) {
 /* ------------------------------------------------------------------ */
 /* decInfinity -- set canonical Infinity with sign from a decFloat    */
 /*								      */
 /*   result gets a canonical Infinity				      */
-/*   df	    is source decFloat (only the sign is used)		      */
+/*   df     is source decFloat (only the sign is used)		      */
 /*   returns result						      */
 /*								      */
-/* df may be the same as 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 */
+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)				      */
+/*   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						      */
 /*								      */
 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		      */
+/* 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		      */
+/*   returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr 	      */
 /*								      */
-/* No error is possible; status and mode are unchanged.		      */
+/* 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 */
+uInt	uiwork; 			/* for macros */
 /* 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;
 /* 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 */
+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 */
 }
 /* decode the coefficients */
 /* (shift both right two if Quad to make a multiple of four) */
 #if QUAD
-ub=bufl;                            /* avoid type-pun violation */
+UBFROMUI(bufl, 0);
-UINTAT(ub)=0;
+UBFROMUI(bufr, 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 */
+uInt ui=UBTOUI(ub);
+if (ui==UBTOUI(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 */
+UBFROMUI(bufl+DECPMAX+QUAD*2, 0);	/* add eight zeros */
-UINTAT(bufl+DECPMAX+QUAD*2+4)=0;	/* .. */
+UBFROMUI(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 *up;			/* work */
 uByte *upend=bufl+DECPMAX+QUAD*2+shift;
-for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
+for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
 /* [pads up to 36 in all for Quad] */
 for (;; ub+=4) {
-	if (UINTAT(ub)!=0) return sigl;
+	if (UBTOUI(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 */
+uInt ui=UBTOUI(ub);
+if (ui!=UBTOUI(uc)) {		/* mismatch found */
 	for (;; uc++, ub++) {		/* check from left [little-endian?] */
 	  if (*ub>*uc) return sigl;	/* difference found */
 	  if (*ub<*uc) return sigr;	/* .. */
 	  }
 	} /* mismatch */
 } /* 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 */
+UBFROMUI(bufr+DECPMAX+QUAD*2, 0);	/* add eight zeros */
-UINTAT(bufr+DECPMAX+QUAD*2+4)=0;	/* .. */
+UBFROMUI(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 *up;			/* work */
 uByte *upend=bufr+DECPMAX+QUAD*2-shift;
-for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UINTAT(up)=0;
+for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
 /* [pads up to 36 in all for Quad] */
 for (;; uc+=4) {
-	if (UINTAT(uc)!=0) return sigr;
+	if (UBTOUI(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 */
+uInt ui=UBTOUI(ub);
+if (ui!=UBTOUI(uc)) {		/* mismatch found */
 	for (;; ub++, uc++) {		/* check from left [little-endian?] */
 	  if (*ub>*uc) return sigl;	/* difference found */
 	  if (*ub<*uc) return sigr;	/* .. */
 	  }
 	} /* mismatch */
 } /* decNumCompare */
 /* ------------------------------------------------------------------ */
 /* decToInt32 -- local routine to effect ToInteger conversions	      */
 /*								      */
-/*   df	    is the decFloat to convert				      */
+/*   df     is the decFloat to convert				      */
 /*   set    is the context					      */
-/*   rmode  is the rounding mode to use				      */
+/*   rmode  is the rounding mode to use 			      */
-/*   exact  is 1 if Inexact should be signalled			      */
+/*   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 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 */
+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;
 }
 /* ------------------------------------------------------------------ */
 /* decToIntegral -- local routine to effect ToIntegral value	      */
 /*								      */
 /*   result gets the result					      */
-/*   df	    is the decFloat to round				      */
+/*   df     is the decFloat to round				      */
 /*   set    is the context					      */
-/*   rmode  is the rounding mode to use				      */
+/*   rmode  is the rounding mode to use 			      */
-/*   exact  is 1 if Inexact should be signalled			      */
+/*   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) {
 enum rounding saveround;	   /* saver */
 uInt savestatus;		   /* .. */
 decFloat zero;		   /* work */
 /* Start decoding the argument */
-sourhi=DFWORD(df, 0);		   /* top word */
+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);
 /* 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 .. */
+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 .. */
+set->round=saveround; 		/* restore rounding mode .. */
 if (!exact) set->status=savestatus;	/* .. and status, unless exact */
 return result;
 } /* decToIntegral */

Mercurial > hg > CbC > CbC_gcc

comparison libdecnumber/decBasic.c @ 55:77e2b8dfacca gcc-4.4.5