/* Common base code for the decNumber C Library.   Copyright (C) 2007-2026 Free Software Foundation, Inc.   Contributed by IBM Corporation.  Author Mike Cowlishaw.    This file is part of GCC.    GCC is free software; you can redistribute it and/or modify it under   the terms of the GNU General Public License as published by the Free   Software Foundation; either version 3, or (at your option) any later   version.    GCC is distributed in the hope that it will be useful, but WITHOUT ANY   WARRANTY; without even the implied warranty of MERCHANTABILITY or   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   for more details. Under Section 7 of GPL version 3, you are granted additionalpermissions described in the GCC Runtime Library Exception, version3.1, as published by the Free Software Foundation. You should have received a copy of the GNU General Public License anda copy of the GCC Runtime Library Exception along with this program;see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see<http://www.gnu.org/licenses/>.  */ /* ------------------------------------------------------------------ *//* decBasic.c -- common base code for Basic decimal types	      *//* ------------------------------------------------------------------ *//* This module comprises code that is shared between decDouble and    *//* decQuad (but not decSingle).  The main arithmetic operations are   *//* here (Add, Subtract, Multiply, FMA, and Division operators).       *//*								      *//* Unlike decNumber, parameterization takes place at compile time     *//* rather than at runtime.  The parameters are set in the decDouble.c *//* (etc.) files, which then include this one to produce the compiled  *//* code.  The functions here, therefore, are code shared between      *//* multiple formats.						      *//*								      *//* This must be included after decCommon.c.			      *//* ------------------------------------------------------------------ *//* Names here refer to decFloat rather than to decDouble, etc., and *//* the functions are in strict alphabetical order. */ /* The compile-time flags SINGLE, DOUBLE, and QUAD are set up in *//* decCommon.c */#if !defined(QUAD)  #error decBasic.c must be included after decCommon.c#endif#if SINGLE  #error Routines in decBasic.c are for decDouble and decQuad only#endif /* Private constants */#define DIVIDE	    0x80000000	   /* Divide operations [as flags] */#define REMAINDER   0x40000000	   /* .. */#define DIVIDEINT   0x20000000	   /* .. */#define REMNEAR     0x10000000	   /* .. */ /* Private functions (local, used only by routines in this module) */static decFloat *decDivide(decFloat *, const decFloat *,			      const decFloat *, decContext *, uInt);static decFloat *decCanonical(decFloat *, const decFloat *);static void	 decFiniteMultiply(bcdnum *, uByte *, const decFloat *,			      const decFloat *);static decFloat *decInfinity(decFloat *, const decFloat *);static decFloat *decInvalid(decFloat *, decContext *);static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *,			      decContext *);static Int	 decNumCompare(const decFloat *, const decFloat *, Flag);static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *,			      enum rounding, Flag);static uInt	 decToInt32(const decFloat *, decContext *, enum rounding,			      Flag, Flag); /* ------------------------------------------------------------------ *//* decCanonical -- copy a decFloat, making canonical		      *//*								      *//*   result gets the canonicalized df				      *//*   df     is the decFloat to copy and make canonical		      *//*   returns result						      *//*								      *//* This is exposed via decFloatCanonical for Double and Quad only.    *//* This works on specials, too; no error or exception is possible.    *//* ------------------------------------------------------------------ */static decFloat * decCanonical(decFloat *result, const decFloat *df) {  uInt encode, precode, dpd;	   /* work */  uInt inword, uoff, canon;	   /* .. */  Int  n;			   /* counter (down) */  if (df!=result) *result=*df;	   /* effect copy if needed */  if (DFISSPECIAL(result)) {    if (DFISINF(result)) return decInfinity(result, df); /* clean Infinity */    /* is a NaN */    DFWORD(result, 0)&=~ECONNANMASK;	/* clear ECON except selector */    if (DFISCCZERO(df)) return result;	/* coefficient continuation is 0 */    /* drop through to check payload */    }  /* return quickly if the coefficient continuation is canonical */  { /* declare block */  #if DOUBLE    uInt sourhi=DFWORD(df, 0);    uInt sourlo=DFWORD(df, 1);    if (CANONDPDOFF(sourhi, 8)     && CANONDPDTWO(sourhi, sourlo, 30)     && CANONDPDOFF(sourlo, 20)     && CANONDPDOFF(sourlo, 10)     && CANONDPDOFF(sourlo, 0)) return result;  #elif QUAD    uInt sourhi=DFWORD(df, 0);    uInt sourmh=DFWORD(df, 1);    uInt sourml=DFWORD(df, 2);    uInt sourlo=DFWORD(df, 3);    if (CANONDPDOFF(sourhi, 4)     && CANONDPDTWO(sourhi, sourmh, 26)     && CANONDPDOFF(sourmh, 16)     && CANONDPDOFF(sourmh, 6)     && CANONDPDTWO(sourmh, sourml, 28)     && CANONDPDOFF(sourml, 18)     && CANONDPDOFF(sourml, 8)     && CANONDPDTWO(sourml, sourlo, 30)     && CANONDPDOFF(sourlo, 20)     && CANONDPDOFF(sourlo, 10)     && CANONDPDOFF(sourlo, 0)) return result;  #endif  } /* block */   /* Loop to repair a non-canonical coefficent, as needed */  inword=DECWORDS-1;		   /* current input word */  uoff=0;			   /* bit offset of declet */  encode=DFWORD(result, inword);  for (n=DECLETS-1; n>=0; n--) {   /* count down declets of 10 bits */    dpd=encode>>uoff;    uoff+=10;    if (uoff>32) {		   /* crossed uInt boundary */      inword--;      encode=DFWORD(result, inword);      uoff-=32;      dpd|=encode<<(10-uoff);	   /* get pending bits */      }    dpd&=0x3ff; 		   /* clear uninteresting bits */    if (dpd<0x16e) continue;	   /* must be canonical */    canon=BIN2DPD[DPD2BIN[dpd]];   /* determine canonical declet */    if (canon==dpd) continue;	   /* have canonical declet */    /* need to replace declet */    if (uoff>=10) {		   /* all within current word */      encode&=~(0x3ff<<(uoff-10)); /* clear the 10 bits ready for replace */      encode|=canon<<(uoff-10);    /* insert the canonical form */      DFWORD(result, inword)=encode;	/* .. and save */      continue;      }    /* straddled words */    precode=DFWORD(result, inword+1);	/* get previous */    precode&=0xffffffff>>(10-uoff);	/* clear top bits */    DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff)));    encode&=0xffffffff<<uoff;		/* clear bottom bits */    encode|=canon>>(10-uoff);		/* insert canonical */    DFWORD(result, inword)=encode;	/* .. and save */    } /* n */  return result;  } /* decCanonical */ /* ------------------------------------------------------------------ *//* decDivide -- divide operations				      *//*								      *//*   result gets the result of dividing dfl by dfr:		      *//*   dfl    is the first decFloat (lhs) 			      *//*   dfr    is the second decFloat (rhs)			      *//*   set    is the context					      *//*   op     is the operation selector				      *//*   returns result						      *//*								      *//* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR.	      *//* ------------------------------------------------------------------ */#define DIVCOUNT  0		   /* 1 to instrument subtractions counter */#define DIVBASE   ((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	 quo[DIVOPLEN+1];	   /* quotient in base-billion .. */  uByte  bcdacc[(DIVOPLEN+1)*9+2]; /* for quotient in BCD, +1, +1 */  uInt	 *msua, *msud, *msuq;	   /* -> msu of acc, div, and quo */  Int	 divunits, accunits;	   /* lengths */  Int	 quodigits;		   /* digits in quotient */  uInt	 *lsua, *lsuq;		   /* -> current acc and quo lsus */  Int	 length, multiplier;	   /* work */  uInt	 carry, sign;		   /* .. */  uInt	 *ua, *ud, *uq; 	   /* .. */  uByte  *ub;			   /* .. */  uInt	 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   /* calculate sign */  num.sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;   if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */    /* NaNs are handled as usual */    if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);    /* one or two infinities */    if (DFISINF(dfl)) {      if (DFISINF(dfr)) return decInvalid(result, set); /* Two infinities bad */      if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* as is rem */      /* Infinity/x is infinite and quiet, even if x=0 */      DFWORD(result, 0)=num.sign;      return decInfinity(result, result);      }    /* must be x/Infinity -- remainders are lhs */    if (op&(REMAINDER|REMNEAR)) return decCanonical(result, dfl);    /* divides: return zero with correct sign and exponent depending */    /* on op (Etiny for divide, 0 for divideInt) */    decFloatZero(result);    if (op==DIVIDEINT) DFWORD(result, 0)|=num.sign; /* add sign */     else DFWORD(result, 0)=num.sign;	     /* zeros the exponent, too */    return result;    }  /* next, handle zero operands (x/0 and 0/x) */  if (DFISZERO(dfr)) {			     /* x/0 */    if (DFISZERO(dfl)) {		     /* 0/0 is undefined */      decFloatZero(result);      DFWORD(result, 0)=DECFLOAT_qNaN;      set->status|=DEC_Division_undefined;      return result;      }    if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* bad rem */    set->status|=DEC_Division_by_zero;    DFWORD(result, 0)=num.sign;    return decInfinity(result, result);      /* x/0 -> signed Infinity */    }  num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr);  /* ideal exponent */  if (DFISZERO(dfl)) {			     /* 0/x (x!=0) */    /* if divide, result is 0 with ideal exponent; divideInt has */    /* exponent=0, remainders give zero with lower exponent */    if (op&DIVIDEINT) {      decFloatZero(result);      DFWORD(result, 0)|=num.sign;	     /* add sign */      return result;      }    if (!(op&DIVIDE)) { 		     /* a remainder */      /* exponent is the minimum of the operands */      num.exponent=MINI(GETEXPUN(dfl), GETEXPUN(dfr));      /* if the result is zero the sign shall be sign of dfl */      num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;      }    bcdacc[0]=0;    num.msd=bcdacc;			     /* -> 0 */    num.lsd=bcdacc;			     /* .. */    return decFinalize(result, &num, set);   /* [divide may clamp exponent] */    } /* 0/x */  /* [here, both operands are known to be finite and non-zero] */   /* extract the operand coefficents into 'units' which are */  /* base-billion; the lhs is high-aligned in acc and the msu of both */  /* acc and div is at the right-hand end of array (offset length-1); */  /* the quotient can need one more unit than the operands as digits */  /* in it are not necessarily aligned neatly; further, the quotient */  /* may not start accumulating until after the end of the initial */  /* operand in acc if that is small (e.g., 1) so the accumulator */  /* must have at least that number of units extra (at the ls end) */  GETCOEFFBILL(dfl, acc+DIVACCLEN-DIVOPLEN);  GETCOEFFBILL(dfr, div);  /* zero the low uInts of acc */  acc[0]=0;  acc[1]=0;  acc[2]=0;  acc[3]=0;  #if DOUBLE    #if DIVOPLEN!=2      #error Unexpected Double DIVOPLEN    #endif  #elif QUAD  acc[4]=0;  acc[5]=0;  acc[6]=0;  acc[7]=0;    #if DIVOPLEN!=4      #error Unexpected Quad DIVOPLEN    #endif  #endif   /* set msu and lsu pointers */  msua=acc+DIVACCLEN-1;       /* [leading zeros removed below] */  msuq=quo+DIVOPLEN;  /*[loop for div will terminate because operands are non-zero] */  for (msud=div+DIVOPLEN-1; *msud==0;) msud--;  /* the initial least-significant unit of acc is set so acc appears */  /* to have the same length as div. */  /* This moves one position towards the least possible for each */  /* iteration */  divunits=(Int)(msud-div+1); /* precalculate */  lsua=msua-divunits+1;       /* initial working lsu of acc */  lsuq=msuq;		      /* and of quo */   /* set up the estimator for the multiplier; this is the msu of div, */  /* plus two bits from the unit below (if any) rounded up by one if */  /* there are any non-zero bits or units below that [the extra two */  /* bits makes for a much better estimate when the top unit is small] */  divtop=*msud<<2;  if (divunits>1) {    uInt *um=msud-1;    uInt d=*um;    if (d>=750000000) {divtop+=3; d-=750000000;}     else if (d>=500000000) {divtop+=2; d-=500000000;}     else if (d>=250000000) {divtop++; d-=250000000;}    if (d) divtop++;     else for (um--; um>=div; um--) if (*um) {      divtop++;      break;      }    } /* >1 unit */   #if DECTRACE  {Int i;  printf("----- div=");  for (i=divunits-1; i>=0; i--) printf("%09ld ", (LI)div[i]);  printf("\n");}  #endif   /* now collect up to DECPMAX+1 digits in the quotient (this may */  /* need OPLEN+1 uInts if unaligned) */  quodigits=0;		      /* no digits yet */  for (;; lsua--) {	      /* outer loop -- each input position */    #if DECCHECK    if (lsua<acc) {      printf("Acc underrun...\n");      break;      }    #endif    #if DECTRACE    printf("Outer: quodigits=%ld acc=", (LI)quodigits);    for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua);    printf("\n");    #endif    *lsuq=0;		      /* default unit result is 0 */    for (;;) {		      /* inner loop -- calculate quotient unit */      /* strip leading zero units from acc (either there initially or */      /* from subtraction below); this may strip all if exactly 0 */      for (; *msua==0 && msua>=lsua;) msua--;      accunits=(Int)(msua-lsua+1);		  /* [maybe 0] */      /* subtraction is only necessary and possible if there are as */      /* least as many units remaining in acc for this iteration as */      /* there are in div */      if (accunits<divunits) {	if (accunits==0) msua++;		  /* restore */	break;	}       /* If acc is longer than div then subtraction is definitely */      /* possible (as msu of both is non-zero), but if they are the */      /* same length a comparison is needed. */      /* If a subtraction is needed then a good estimate of the */      /* multiplier for the subtraction is also needed in order to */      /* minimise the iterations of this inner loop because the */      /* subtractions needed dominate division performance. */      if (accunits==divunits) {	/* compare the high divunits of acc and div: */	/* acc<div:  this quotient unit is unchanged; subtraction */	/*	     will be possible on the next iteration */	/* acc==div: quotient gains 1, set acc=0 */	/* acc>div:  subtraction necessary at this position */	for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break;	/* [now at first mismatch or lsu] */	if (*ud>*ua) break;			  /* next time... */	if (*ud==*ua) { 			  /* all compared equal */	  *lsuq+=1;				  /* increment result */	  msua=lsua;				  /* collapse acc units */	  *msua=0;				  /* .. to a zero */	  break;	  } 	/* subtraction necessary; estimate multiplier [see above] */	/* if both *msud and *msua are small it is cost-effective to */	/* bring in part of the following units (if any) to get a */	/* better estimate (assume some other non-zero in div) */	#define DIVLO 1000000U	#define DIVHI (DIVBASE/DIVLO)	#if DECUSE64	  if (divunits>1) {	    /* there cannot be a *(msud-2) for DECDOUBLE so next is */	    /* an exact calculation unless DECQUAD (which needs to */	    /* assume bits out there if divunits>2) */	    uLong mul=(uLong)*msua * DIVBASE + *(msua-1);	    uLong div=(uLong)*msud * DIVBASE + *(msud-1);	    #if QUAD	    if (divunits>2) div++;	    #endif	    mul/=div;	    multiplier=(Int)mul;	    }	   else multiplier=*msua/(*msud);	#else	  if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {	    multiplier=(*msua*DIVHI + *(msua-1)/DIVLO)		      /(*msud*DIVHI + *(msud-1)/DIVLO +1);	    }	   else multiplier=(*msua<<2)/divtop;	#endif	}       else {					  /* accunits>divunits */	/* msud is one unit 'lower' than msua, so estimate differently */	#if DECUSE64	  uLong mul;	  /* as before, bring in extra digits if possible */	  if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {	    mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI	       + *(msua-2)/DIVLO;	    mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1);	    }	   else if (divunits==1) {	    mul=(uLong)*msua * DIVBASE + *(msua-1);	    mul/=*msud;       /* no more to the right */	    }	   else {	    mul=(uLong)(*msua) * (uInt)(DIVBASE<<2)		+ (*(msua-1)<<2);	    mul/=divtop;      /* [divtop already allows for sticky bits] */	    }	  multiplier=(Int)mul;	#else	  multiplier=*msua * ((DIVBASE<<2)/divtop);	#endif	}      if (multiplier==0) multiplier=1;		  /* marginal case */      *lsuq+=multiplier;       #if DIVCOUNT      /* printf("Multiplier: %ld\n", (LI)multiplier); */      divcount++;      #endif       /* Carry out the subtraction  acc-(div*multiplier); for each */      /* unit in div, do the multiply, split to units (see */      /* decFloatMultiply for the algorithm), and subtract from acc */      #define DIVMAGIC	2305843009U		  /* 2**61/10**9 */      #define DIVSHIFTA 29      #define DIVSHIFTB 32      carry=0;      for (ud=div, ua=lsua; ud<=msud; ud++, ua++) {	uInt lo, hop;	#if DECUSE64	  uLong sub=(uLong)multiplier*(*ud)+carry;	  if (sub<DIVBASE) {	    carry=0;	    lo=(uInt)sub;	    }	   else {	    hop=(uInt)(sub>>DIVSHIFTA);	    carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB);	    /* the estimate is now in hi; now calculate sub-hi*10**9 */	    /* to get the remainder (which will be <DIVBASE)) */	    lo=(uInt)sub;	    lo-=carry*DIVBASE;			  /* low word of result */	    if (lo>=DIVBASE) {	      lo-=DIVBASE;			  /* correct by +1 */	      carry++;	      }	    }	#else /* 32-bit */	  uInt hi;	  /* calculate multiplier*(*ud) into hi and lo */	  LONGMUL32HI(hi, *ud, multiplier);	  /* get the high word */	  lo=multiplier*(*ud);			  /* .. and the low */	  lo+=carry;				  /* add the old hi */	  carry=hi+(lo<carry);			  /* .. with any carry */	  if (carry || lo>=DIVBASE) {		  /* split is needed */	    hop=(carry<<3)+(lo>>DIVSHIFTA);	  /* hi:lo/2**29 */	    LONGMUL32HI(carry, hop, DIVMAGIC);	  /* only need the high word */	    /* [DIVSHIFTB is 32, so carry can be used directly] */	    /* the estimate is now in carry; now calculate hi:lo-est*10**9; */	    /* happily the top word of the result is irrelevant because it */	    /* will always be zero so this needs only one multiplication */	    lo-=(carry*DIVBASE);	    /* the correction here will be at most +1; do it */	    if (lo>=DIVBASE) {	      lo-=DIVBASE;	      carry++;	      }	    }	#endif	if (lo>*ua) {		   /* borrow needed */	  *ua+=DIVBASE;	  carry++;	  }	*ua-=lo;	} /* ud loop */      if (carry) *ua-=carry;	   /* accdigits>divdigits [cannot borrow] */      } /* inner loop */     /* the outer loop terminates when there is either an exact result */    /* or enough digits; first update the quotient digit count and */    /* pointer (if any significant digits) */    #if DECTRACE    if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq);    #endif

Lines 1–500 of 3909. Load more