Lines Matching +full:remove +full:- +full:item
29 02110-1301, USA. */
31 /* ------------------------------------------------------------------ */
33 /* ------------------------------------------------------------------ */
44 /* If DECDPUN>4 or DECUSE64=1, the C99 64-bit int64_t and */
54 /* range -999,999,999 through 0). Mathematical functions (for */
56 /* tightly: digits, emax, and -emin in the context must be <= */
71 /* SIGFPE signal is then raised if the corresponding trap-enabler */
82 /* representation as it comprises fields which may be machine- */
89 /* Results are undefined if a badly-formed structure (or a NULL */
95 /* given well-formed operands, even if the value of the operands */
100 /* ------------------------------------------------------------------ */
107 /* code which is protected -- break may be used to exit this. */
116 /* action in a routine, as non-0 status may raise a trap and hence */
123 /* 'top level' routines which could cause a non-returning */
132 /* require more careful checks because of the risk of 31-bit */
133 /* overflow (the most negative valid exponent is -1999999997, for */
134 /* a 999999999-digit number with adjusted exponent of -999999999). */
138 /* residue (in the range -1 through 9). This avoids any double- */
149 /* 7. The multiply-by-reciprocal 'trick' is used for partitioning */
153 /* a divide (unless a floating-point or 64-bit multiply is */
157 /* lhs -- left hand side (operand, of an operation) */
158 /* lsd -- least significant digit (of coefficient) */
159 /* lsu -- least significant Unit (of coefficient) */
160 /* msd -- most significant digit (of coefficient) */
161 /* msi -- most significant item (in an array) */
162 /* msu -- most significant Unit (of coefficient) */
163 /* rhs -- right hand side (operand, of an operation) */
164 /* +ve -- positive */
165 /* -ve -- negative */
166 /* ** -- raise to the power */
167 /* ------------------------------------------------------------------ */
170 #include "qemu/host-utils.h"
197 #define BADINT (Int)0x80000000 /* most-negative Int; error indicator */
204 /* Granularity-dependent code */
208 /* Constant multipliers for divide-by-power-of five using reciprocal */
212 /* QUOT10 -- macro to return the quotient of unit u divided by 10**n */
215 /* For DECDPUN>4 non-ANSI-89 64-bit types are needed. */
278 /* masked special-values bits */
279 #define SPECIALARG (rhs->bits & DECSPECIAL)
280 #define SPECIALARGS ((lhs->bits | rhs->bits) & DECSPECIAL)
298 /* if the testing is done in a multi-thread environment. */
308 /* complement of digits (where appropriate -- this is not the case */
329 /* ------------------------------------------------------------------ */
330 /* from-int32 -- conversion from Int or uInt */
337 /* ------------------------------------------------------------------ */
343 else unsig=-in; /* invert */ in decNumberFromInt32()
347 if (in<0) dn->bits=DECNEG; /* sign needed */ in decNumberFromInt32()
355 for (up=dn->lsu; uin>0; up++) { in decNumberFromUInt32()
359 dn->digits=decGetDigits(dn->lsu, up-dn->lsu); in decNumberFromUInt32()
363 /* ------------------------------------------------------------------ */
364 /* to-int32 -- conversion to Int or uInt */
371 /* it is a NaN, Infinite, or out-of-range. */
372 /* ------------------------------------------------------------------ */
379 if (dn->bits&DECSPECIAL || dn->digits>10 || dn->exponent!=0) ; /* bad */ in decNumberToInt32()
384 up=dn->lsu; /* -> lsu */ in decNumberToInt32()
392 for (d=DECDPUN; d<dn->digits; up++, d+=DECDPUN) hi+=*up*powers[d-1]; in decNumberToInt32()
395 /* most-negative is a reprieve */ in decNumberToInt32()
396 if (dn->bits&DECNEG && hi==214748364 && lo==8) return 0x80000000; in decNumberToInt32()
397 /* bad -- drop through */ in decNumberToInt32()
399 else { /* in-range always */ in decNumberToInt32()
401 if (dn->bits&DECNEG) return -i; in decNumberToInt32()
414 if (dn->bits&DECSPECIAL || dn->digits>10 || dn->exponent!=0 in decNumberToUInt32()
415 || (dn->bits&DECNEG && !ISZERO(dn))); /* bad */ in decNumberToUInt32()
420 up=dn->lsu; /* -> lsu */ in decNumberToUInt32()
428 for (d=DECDPUN; d<dn->digits; up++, d+=DECDPUN) hi+=*up*powers[d-1]; in decNumberToUInt32()
442 unsig = -unsig; in decNumberFromInt64()
447 dn->bits = DECNEG; /* sign needed */ in decNumberFromInt64()
459 for (up = dn->lsu; uin > 0; up++) { in decNumberFromUInt64()
463 dn->digits = decGetDigits(dn->lsu, up-dn->lsu); in decNumberFromUInt64()
472 unsig_hi = -unsig_hi; in decNumberFromInt128()
475 lo = -lo; in decNumberFromInt128()
481 dn->bits = DECNEG; /* sign needed */ in decNumberFromInt128()
494 for (up = dn->lsu; hi > 0 || lo > 0; up++) { in decNumberFromUInt128()
498 dn->digits = decGetDigits(dn->lsu, up - dn->lsu); in decNumberFromUInt128()
502 /* ------------------------------------------------------------------ */
503 /* to-int64 -- conversion to int64 */
512 /* ------------------------------------------------------------------ */
516 if (decNumberIsSpecial(dn) || (dn->exponent < 0) || in decNumberIntegralToInt64()
517 (dn->digits + dn->exponent > 19)) { in decNumberIntegralToInt64()
523 up = dn->lsu; /* -> lsu */ in decNumberIntegralToInt64()
525 for (d = 1; d <= dn->digits; up++, d += DECDPUN) { in decNumberIntegralToInt64()
527 hi += *up * powers[d-1]; in decNumberIntegralToInt64()
534 hi *= (uint64_t)powers[dn->exponent]; in decNumberIntegralToInt64()
538 return (decNumberIsNegative(dn)) ? -((int64_t)hi) : (int64_t)hi; in decNumberIntegralToInt64()
546 /* ------------------------------------------------------------------ */
547 /* decNumberIntegralToInt128 -- conversion to int128 */
556 /* ------------------------------------------------------------------ */
565 if (decNumberIsSpecial(dn) || (dn->exponent < 0) || in decNumberIntegralToInt128()
566 (dn->digits + dn->exponent > 39)) { in decNumberIntegralToInt128()
570 up = dn->lsu; /* -> lsu */ in decNumberIntegralToInt128()
572 for (d = (dn->digits - 1) / DECDPUN; d >= 0; d--) { in decNumberIntegralToInt128()
585 if (mulUInt128ByPowOf10(&lo, &hi, dn->exponent)) { in decNumberIntegralToInt128()
592 *phigh = -hi; in decNumberIntegralToInt128()
596 *plow = -lo; in decNumberIntegralToInt128()
609 /* ------------------------------------------------------------------ */
610 /* to-scientific-string -- conversion to numeric string */
611 /* to-engineering-string -- conversion to numeric string */
619 /* string must be at least dn->digits+14 characters long */
622 /* ------------------------------------------------------------------ */
633 /* ------------------------------------------------------------------ */
634 /* to-number -- conversion from numeric string */
636 /* decNumberFromString -- convert string to decNumber */
637 /* dn -- the number structure to fill */
638 /* chars[] -- the string to convert ('\0' terminated) */
639 /* set -- the context used for processing any error, */
652 /* ------------------------------------------------------------------ */
660 Unit *allocres=NULL; /* -> allocated result, iff allocated */ in decNumberFromString()
663 const char *cfirst=chars; /* -> first character of decimal part */ in decNumberFromString()
664 const char *last=NULL; /* -> last digit of decimal part */ in decNumberFromString()
679 for (c=chars;; c++) { /* -> input character */ in decNumberFromString()
690 if (*c=='-') { /* valid - sign */ in decNumberFromString()
698 /* *c is not a digit, or a valid +, -, or '.' */ in decNumberFromString()
707 if (!set->extended) break; /* hopeless */ in decNumberFromString()
714 dn->bits=bits | DECINF; in decNumberFromString()
720 dn->bits=bits | DECNAN; /* assume simple NaN */ in decNumberFromString()
723 dn->bits=bits | DECSNAN; in decNumberFromString()
732 /* -> start of integer and skip leading 0s [including plain 0] */ in decNumberFromString()
744 if (d>set->digits-1) { in decNumberFromString()
746 /* clamped, in which case can only be digits-1] */ in decNumberFromString()
747 if (set->clamp) break; in decNumberFromString()
748 if (d>set->digits) break; in decNumberFromString()
752 bits=dn->bits; /* for copy-back */ in decNumberFromString()
758 const char *firstexp; /* -> first significant exponent digit */ in decNumberFromString()
761 /* Found 'e' or 'E' -- now process explicit exponent */ in decNumberFromString()
765 if (*c=='-') {nege=1; c++;} in decNumberFromString()
773 exponent=X10(exponent)+(Int)*c-(Int)'0'; in decNumberFromString()
783 /* [up to 1999999999 is OK, for example 1E-1000000998] */ in decNumberFromString()
785 if (nege) exponent=-exponent; /* was negative */ in decNumberFromString()
790 /* cfirst->first digit (never dot), last->last digit (ditto) */ in decNumberFromString()
796 if (*c!='0') break; /* non-zero found */ in decNumberFromString()
797 d--; /* 0 stripped */ in decNumberFromString()
801 if (*cfirst=='0' && !set->extended) { in decNumberFromString()
809 if (dotchar!=NULL && dotchar<last) /* non-trailing '.' found? */ in decNumberFromString()
810 exponent-=(last-dotchar); /* adjust exponent */ in decNumberFromString()
815 if (d<=set->digits) res=dn->lsu; /* fits into supplied decNumber */ in decNumberFromString()
825 /* res now -> number lsu, buffer, or allocated storage for Unit array */ in decNumberFromString()
831 up=res+D2U(d)-1; /* -> msu */ in decNumberFromString()
832 cut=d-(up-res)*DECDPUN; /* digits in top unit */ in decNumberFromString()
835 out=X10(out)+(Int)*c-(Int)'0'; in decNumberFromString()
837 cut--; in decNumberFromString()
840 up--; /* prepare for unit below.. */ in decNumberFromString()
848 up=res; /* -> lsu */ in decNumberFromString()
849 for (c=last; c>=cfirst; c--) { /* over each character, from least */ in decNumberFromString()
851 *up=(Unit)((Int)*c-(Int)'0'); in decNumberFromString()
856 dn->bits=bits; in decNumberFromString()
857 dn->exponent=exponent; in decNumberFromString()
858 dn->digits=d; in decNumberFromString()
861 if (d>set->digits) { in decNumberFromString()
869 if ((dn->exponent-1<set->emin-dn->digits) in decNumberFromString()
870 || (dn->exponent-1>set->emax-set->digits)) { in decNumberFromString()
887 /* ------------------------------------------------------------------ */
888 /* decNumberAbs -- absolute value operator */
897 /* C must have space for set->digits digits. */
898 /* ------------------------------------------------------------------ */
901 /* ------------------------------------------------------------------ */
912 dzero.exponent=rhs->exponent; /* [no coefficient expansion] */ in decNumberAbs()
913 decAddOp(res, &dzero, rhs, set, (uByte)(rhs->bits & DECNEG), &status); in decNumberAbs()
921 /* ------------------------------------------------------------------ */
922 /* decNumberAdd -- add two Numbers */
931 /* C must have space for set->digits digits. */
932 /* ------------------------------------------------------------------ */
945 /* ------------------------------------------------------------------ */
946 /* decNumberAnd -- AND two Numbers, digitwise */
955 /* C must have space for set->digits digits. */
959 /* ------------------------------------------------------------------ */
962 const Unit *ua, *ub; /* -> operands */ in decNumberAnd()
963 const Unit *msua, *msub; /* -> operand msus */ in decNumberAnd()
964 Unit *uc, *msuc; /* -> result and its msu */ in decNumberAnd()
970 if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs) in decNumberAnd()
971 || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) { in decNumberAnd()
977 ua=lhs->lsu; /* bottom-up */ in decNumberAnd()
978 ub=rhs->lsu; /* .. */ in decNumberAnd()
979 uc=res->lsu; /* .. */ in decNumberAnd()
980 msua=ua+D2U(lhs->digits)-1; /* -> msu of lhs */ in decNumberAnd()
981 msub=ub+D2U(rhs->digits)-1; /* -> msu of rhs */ in decNumberAnd()
982 msuc=uc+D2U(set->digits)-1; /* -> msu of result */ in decNumberAnd()
983 msudigs=MSUDIGITS(set->digits); /* [faster than remainder] */ in decNumberAnd()
1005 if (uc==msuc && i==msudigs-1) break; /* just did final digit */ in decNumberAnd()
1009 /* [here uc-1 is the msu of the result] */ in decNumberAnd()
1010 res->digits=decGetDigits(res->lsu, uc-res->lsu); in decNumberAnd()
1011 res->exponent=0; /* integer */ in decNumberAnd()
1012 res->bits=0; /* sign=0 */ in decNumberAnd()
1016 /* ------------------------------------------------------------------ */
1017 /* decNumberCompare -- compare two Numbers */
1027 /* ------------------------------------------------------------------ */
1036 /* ------------------------------------------------------------------ */
1037 /* decNumberCompareSignal -- compare, signalling on all NaNs */
1047 /* ------------------------------------------------------------------ */
1056 /* ------------------------------------------------------------------ */
1057 /* decNumberCompareTotal -- compare two Numbers, using total ordering */
1067 /* -1, 0, or 1. */
1068 /* ------------------------------------------------------------------ */
1077 /* ------------------------------------------------------------------ */
1078 /* decNumberCompareTotalMag -- compare, total ordering of magnitudes */
1088 /* -1, 0, or 1. */
1089 /* ------------------------------------------------------------------ */
1095 decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */ in decNumberCompareTotalMag()
1097 decNumber *allocbufb=NULL; /* -> allocated bufb, iff allocated */ in decNumberCompareTotalMag()
1108 needbytes=sizeof(decNumber)+(D2U(lhs->digits)-1)*sizeof(Unit); in decNumberCompareTotalMag()
1111 if (allocbufa==NULL) { /* hopeless -- abandon */ in decNumberCompareTotalMag()
1117 a->bits&=~DECNEG; /* .. and clear the sign */ in decNumberCompareTotalMag()
1122 needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit); in decNumberCompareTotalMag()
1125 if (allocbufb==NULL) { /* hopeless -- abandon */ in decNumberCompareTotalMag()
1131 b->bits&=~DECNEG; /* .. and clear the sign */ in decNumberCompareTotalMag()
1143 /* ------------------------------------------------------------------ */
1144 /* decNumberDivide -- divide one number by another */
1153 /* C must have space for set->digits digits. */
1154 /* ------------------------------------------------------------------ */
1166 /* ------------------------------------------------------------------ */
1167 /* decNumberDivideInteger -- divide and return integer quotient */
1176 /* C must have space for set->digits digits. */
1177 /* ------------------------------------------------------------------ */
1186 /* ------------------------------------------------------------------ */
1187 /* decNumberExp -- exponentiation */
1195 /* C must have space for set->digits digits. */
1201 /* when A is a zero or -Infinity (giving 1 or 0 respectively). */
1206 /* ------------------------------------------------------------------ */
1209 /* exp(-a) where a can be the tiniest number (Ntiny). */
1210 /* ------------------------------------------------------------------ */
1215 decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */ in decNumberExp()
1228 if (!set->extended) { in decNumberExp()
1230 if (rhs->digits>set->digits) { in decNumberExp()
1251 /* ------------------------------------------------------------------ */
1252 /* decNumberFMA -- fused multiply add */
1265 /* C must have space for set->digits digits. */
1266 /* ------------------------------------------------------------------ */
1274 decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */ in decNumberFMA()
1285 if (!set->extended) { /* [undefined if subset] */ in decNumberFMA()
1295 dcmul.digits=lhs->digits+rhs->digits; /* just enough */ in decNumberFMA()
1296 /* [The above may be an over-estimate for subset arithmetic, but that's OK] */ in decNumberFMA()
1301 needbytes=sizeof(decNumber)+(D2U(dcmul.digits)-1)*sizeof(Unit); in decNumberFMA()
1304 if (allocbufa==NULL) { /* hopeless -- abandon */ in decNumberFMA()
1321 res->bits=DECNAN; in decNumberFMA()
1324 decNumberZero(&dzero); /* make 0 (any non-NaN would do) */ in decNumberFMA()
1332 /* add the third operand and result -> res, and all is done */ in decNumberFMA()
1344 /* ------------------------------------------------------------------ */
1345 /* decNumberInvert -- invert a Number, digitwise */
1353 /* C must have space for set->digits digits. */
1357 /* ------------------------------------------------------------------ */
1360 const Unit *ua, *msua; /* -> operand and its msu */ in decNumberInvert()
1361 Unit *uc, *msuc; /* -> result and its msu */ in decNumberInvert()
1367 if (rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) { in decNumberInvert()
1372 ua=rhs->lsu; /* bottom-up */ in decNumberInvert()
1373 uc=res->lsu; /* .. */ in decNumberInvert()
1374 msua=ua+D2U(rhs->digits)-1; /* -> msu of rhs */ in decNumberInvert()
1375 msuc=uc+D2U(set->digits)-1; /* -> msu of result */ in decNumberInvert()
1376 msudigs=MSUDIGITS(set->digits); /* [faster than remainder] */ in decNumberInvert()
1393 if (uc==msuc && i==msudigs-1) break; /* just did final digit */ in decNumberInvert()
1396 /* [here uc-1 is the msu of the result] */ in decNumberInvert()
1397 res->digits=decGetDigits(res->lsu, uc-res->lsu); in decNumberInvert()
1398 res->exponent=0; /* integer */ in decNumberInvert()
1399 res->bits=0; /* sign=0 */ in decNumberInvert()
1403 /* ------------------------------------------------------------------ */
1404 /* decNumberLn -- natural logarithm */
1412 /* C must have space for set->digits digits. */
1415 /* A<0 -> Invalid */
1416 /* A=0 -> -Infinity (Exact) */
1417 /* A=+Infinity -> +Infinity (Exact) */
1418 /* A=1 exactly -> 0 (Exact) */
1426 /* ------------------------------------------------------------------ */
1430 /* ------------------------------------------------------------------ */
1435 decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */ in decNumberLn()
1446 if (!set->extended) { in decNumberLn()
1448 if (rhs->digits>set->digits) { in decNumberLn()
1454 if (ISZERO(rhs)) { /* +/- zeros -> error */ in decNumberLn()
1473 /* ------------------------------------------------------------------ */
1474 /* decNumberLogB - get adjusted exponent, by 754r rules */
1484 /* -1999999999). */
1487 /* with zeros on the right to set->digits digits while keeping the */
1491 /* A<0 -> Use |A| */
1492 /* A=0 -> -Infinity (Division by zero) */
1493 /* A=Infinite -> +Infinity (Exact) */
1494 /* A=1 exactly -> 0 (Exact) */
1496 /* ------------------------------------------------------------------ */
1505 /* NaNs as usual; Infinities return +Infinity; 0->oops */ in decNumberLogB()
1510 res->bits=DECNEG|DECINF; /* -Infinity */ in decNumberLogB()
1513 else { /* finite non-zero */ in decNumberLogB()
1514 Int ae=rhs->exponent+rhs->digits-1; /* adjusted exponent */ in decNumberLogB()
1522 /* ------------------------------------------------------------------ */
1523 /* decNumberLog10 -- logarithm in base 10 */
1531 /* C must have space for set->digits digits. */
1534 /* A<0 -> Invalid */
1535 /* A=0 -> -Infinity (Exact) */
1536 /* A=+Infinity -> +Infinity (Exact) */
1537 /* A=10**n (if n is an integer) -> n (Exact) */
1545 /* ------------------------------------------------------------------ */
1547 /* ln(A) this is the max(p, rhs->digits + t) + 3, where p is the */
1552 /* ------------------------------------------------------------------ */
1563 decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */ in decNumberLog10()
1566 decNumber *allocbufb=NULL; /* -> allocated bufb, iff allocated */ in decNumberLog10()
1568 decNumber bufw[D2N(10)]; /* working 2-10 digit number */ in decNumberLog10()
1571 decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */ in decNumberLog10()
1584 if (!set->extended) { in decNumberLog10()
1586 if (rhs->digits>set->digits) { in decNumberLog10()
1592 if (ISZERO(rhs)) { /* +/- zeros -> error */ in decNumberLog10()
1601 if (!(rhs->bits&(DECNEG|DECSPECIAL)) && !ISZERO(rhs)) { in decNumberLog10()
1609 if (!(copystat&DEC_Inexact) && w->lsu[0]==1) { in decNumberLog10()
1614 decNumberFromInt32(w, w->exponent); in decNumberLog10()
1622 /* simplify the information-content calculation to use 'total */ in decNumberLog10()
1629 p=(rhs->digits+t>set->digits?rhs->digits+t:set->digits)+3; in decNumberLog10()
1630 needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit); in decNumberLog10()
1633 if (allocbufa==NULL) { /* hopeless -- abandon */ in decNumberLog10()
1640 aset.emin=-DEC_MAX_MATH; /* .. */ in decNumberLog10()
1647 if (a->bits&DECSPECIAL || ISZERO(a)) { in decNumberLog10()
1652 p=set->digits+3; in decNumberLog10()
1653 needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit); in decNumberLog10()
1656 if (allocbufb==NULL) { /* hopeless -- abandon */ in decNumberLog10()
1663 w->lsu[1]=1; w->lsu[0]=0; /* .. */ in decNumberLog10()
1665 w->lsu[0]=10; /* .. */ in decNumberLog10()
1667 w->digits=2; /* .. */ in decNumberLog10()
1672 aset.digits=set->digits; /* for final divide */ in decNumberLog10()
1689 /* ------------------------------------------------------------------ */
1690 /* decNumberMax -- compare two Numbers and return the maximum */
1699 /* C must have space for set->digits digits. */
1700 /* ------------------------------------------------------------------ */
1712 /* ------------------------------------------------------------------ */
1713 /* decNumberMaxMag -- compare and return the maximum by magnitude */
1722 /* C must have space for set->digits digits. */
1723 /* ------------------------------------------------------------------ */
1735 /* ------------------------------------------------------------------ */
1736 /* decNumberMin -- compare two Numbers and return the minimum */
1745 /* C must have space for set->digits digits. */
1746 /* ------------------------------------------------------------------ */
1758 /* ------------------------------------------------------------------ */
1759 /* decNumberMinMag -- compare and return the minimum by magnitude */
1768 /* C must have space for set->digits digits. */
1769 /* ------------------------------------------------------------------ */
1781 /* ------------------------------------------------------------------ */
1782 /* decNumberMinus -- prefix minus operator */
1784 /* This computes C = 0 - A */
1791 /* C must have space for set->digits digits. */
1792 /* ------------------------------------------------------------------ */
1794 /* ------------------------------------------------------------------ */
1805 dzero.exponent=rhs->exponent; /* [no coefficient expansion] */ in decNumberMinus()
1814 /* ------------------------------------------------------------------ */
1815 /* decNumberNextMinus -- next towards -Infinity */
1817 /* This computes C = A - infinitesimal, rounded towards -Infinity */
1824 /* ------------------------------------------------------------------ */
1835 if ((rhs->bits&(DECINF|DECNEG))==DECINF) { in decNumberNextMinus()
1842 dtiny.exponent=DEC_MIN_EMIN-1; /* .. smaller than tiniest */ in decNumberNextMinus()
1850 /* ------------------------------------------------------------------ */
1851 /* decNumberNextPlus -- next towards +Infinity */
1860 /* ------------------------------------------------------------------ */
1870 /* -Infinity is the special case */ in decNumberNextPlus()
1871 if ((rhs->bits&(DECINF|DECNEG))==(DECINF|DECNEG)) { in decNumberNextPlus()
1873 res->bits=DECNEG; /* negative */ in decNumberNextPlus()
1879 dtiny.exponent=DEC_MIN_EMIN-1; /* .. smaller than tiniest */ in decNumberNextPlus()
1887 /* ------------------------------------------------------------------ */
1888 /* decNumberNextToward -- next towards rhs */
1890 /* This computes C = A +/- infinitesimal, rounded towards */
1891 /* +/-Infinity in the direction of B, as per 754r nextafter rules */
1899 /* ------------------------------------------------------------------ */
1921 /* -Infinity is the special case */ in decNumberNextToward()
1922 if ((lhs->bits&(DECINF|DECNEG))==(DECINF|DECNEG)) { in decNumberNextToward()
1924 res->bits=DECNEG; /* negative */ in decNumberNextToward()
1932 if ((lhs->bits&(DECINF|DECNEG))==DECINF) { in decNumberNextToward()
1941 dtiny.exponent=DEC_MIN_EMIN-1; /* .. smaller than tiniest */ in decNumberNextToward()
1942 decAddOp(res, lhs, &dtiny, &workset, sub, &status); /* + or - */ in decNumberNextToward()
1953 /* ------------------------------------------------------------------ */
1954 /* decNumberOr -- OR two Numbers, digitwise */
1963 /* C must have space for set->digits digits. */
1967 /* ------------------------------------------------------------------ */
1970 const Unit *ua, *ub; /* -> operands */ in decNumberOr()
1971 const Unit *msua, *msub; /* -> operand msus */ in decNumberOr()
1972 Unit *uc, *msuc; /* -> result and its msu */ in decNumberOr()
1978 if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs) in decNumberOr()
1979 || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) { in decNumberOr()
1984 ua=lhs->lsu; /* bottom-up */ in decNumberOr()
1985 ub=rhs->lsu; /* .. */ in decNumberOr()
1986 uc=res->lsu; /* .. */ in decNumberOr()
1987 msua=ua+D2U(lhs->digits)-1; /* -> msu of lhs */ in decNumberOr()
1988 msub=ub+D2U(rhs->digits)-1; /* -> msu of rhs */ in decNumberOr()
1989 msuc=uc+D2U(set->digits)-1; /* -> msu of result */ in decNumberOr()
1990 msudigs=MSUDIGITS(set->digits); /* [faster than remainder] */ in decNumberOr()
2011 if (uc==msuc && i==msudigs-1) break; /* just did final digit */ in decNumberOr()
2013 } /* non-zero */ in decNumberOr()
2015 /* [here uc-1 is the msu of the result] */ in decNumberOr()
2016 res->digits=decGetDigits(res->lsu, uc-res->lsu); in decNumberOr()
2017 res->exponent=0; /* integer */ in decNumberOr()
2018 res->bits=0; /* sign=0 */ in decNumberOr()
2022 /* ------------------------------------------------------------------ */
2023 /* decNumberPlus -- prefix plus operator */
2032 /* C must have space for set->digits digits. */
2033 /* ------------------------------------------------------------------ */
2037 /* ------------------------------------------------------------------ */
2047 dzero.exponent=rhs->exponent; /* [no coefficient expansion] */ in decNumberPlus()
2056 /* ------------------------------------------------------------------ */
2057 /* decNumberMultiply -- multiply two Numbers */
2066 /* C must have space for set->digits digits. */
2067 /* ------------------------------------------------------------------ */
2079 /* ------------------------------------------------------------------ */
2080 /* decNumberPower -- raise a number to a power */
2089 /* C must have space for set->digits digits. */
2104 /* ------------------------------------------------------------------ */
2108 decNumber *alloclhs=NULL; /* non-NULL if rounded lhs allocated */ in decNumberPower()
2111 decNumber *allocdac=NULL; /* -> allocated acc buffer, iff used */ in decNumberPower()
2112 decNumber *allocinv=NULL; /* -> allocated 1/x buffer, iff used */ in decNumberPower()
2113 Int reqdigits=set->digits; /* requested DIGITS */ in decNumberPower()
2131 decNumber *dac=dacbuff; /* -> result accumulator */ in decNumberPower()
2141 if (!set->extended) { /* reduce operands and set status, as needed */ in decNumberPower()
2142 if (lhs->digits>reqdigits) { in decNumberPower()
2147 if (rhs->digits>reqdigits) { in decNumberPower()
2162 Flag rhsneg=rhs->bits&DECNEG; /* save rhs sign */ in decNumberPower()
2172 if (rhsneg) res->bits|=DECINF; /* +Infinity [else is +0] */ in decNumberPower()
2174 else if (dac->lsu[0]==0) { /* lhs=1 */ in decNumberPower()
2175 /* 1**Infinity is inexact, so return fully-padded 1.0000 */ in decNumberPower()
2176 Int shift=set->digits-1; in decNumberPower()
2177 *res->lsu=1; /* was 0, make int 1 */ in decNumberPower()
2178 res->digits=decShiftToMost(res->lsu, 1, shift); in decNumberPower()
2179 res->exponent=-shift; /* make 1.0000... */ in decNumberPower()
2183 if (!rhsneg) res->bits|=DECINF; /* +Infinity [else is +0] */ in decNumberPower()
2199 if (decNumberIsNegative(lhs) /* -x .. */ in decNumberPower()
2204 uByte rbits=rhs->bits; /* save */ in decNumberPower()
2206 if (n==0) *res->lsu=1; /* [-]Inf**0 => 1 */ in decNumberPower()
2208 /* -Inf**nonint -> error */ in decNumberPower()
2210 status|=DEC_Invalid_operation; /* -Inf**nonint is error */ in decNumberPower()
2212 if (!(rbits & DECNEG)) bits|=DECINF; /* was not a **-n */ in decNumberPower()
2213 /* [otherwise will be 0 or -0] */ in decNumberPower()
2214 res->bits=bits; in decNumberPower()
2222 if (!set->extended) { /* [unless subset] */ in decNumberPower()
2224 *res->lsu=1; /* return 1 */ in decNumberPower()
2230 uByte rbits=rhs->bits; /* save */ in decNumberPower()
2231 if (rbits & DECNEG) { /* was a 0**(-n) */ in decNumberPower()
2233 if (!set->extended) { /* [bad if subset] */ in decNumberPower()
2240 /* [otherwise will be 0 or -0] */ in decNumberPower()
2241 res->bits=bits; in decNumberPower()
2246 /* integer path. Next handle the non-integer cases */ in decNumberPower()
2247 if (!useint) { /* non-integral rhs */ in decNumberPower()
2248 /* any -ve lhs is bad, as is either operand or context out of */ in decNumberPower()
2258 aset.emin=-DEC_MAX_MATH; /* .. */ in decNumberPower()
2271 aset.digits=MAXI(lhs->digits, set->digits)+6+4; in decNumberPower()
2272 } /* non-integer rhs */ in decNumberPower()
2274 else { /* rhs is in-range integer */ in decNumberPower()
2278 *res->lsu=1; /* .. */ in decNumberPower()
2280 /* rhs is a non-zero integer */ in decNumberPower()
2281 if (n<0) n=-n; /* use abs(n) */ in decNumberPower()
2286 aset.digits=reqdigits+(rhs->digits+rhs->exponent)+2; in decNumberPower()
2288 if (!set->extended) aset.digits--; /* use classic precision */ in decNumberPower()
2296 needbytes=sizeof(decNumber)+(D2U(aset.digits)-1)*sizeof(Unit); in decNumberPower()
2300 if (allocdac==NULL) { /* hopeless -- abandon */ in decNumberPower()
2307 if (!useint) { /* non-integral rhs */ in decNumberPower()
2308 /* x ** y; special-case x=1 here as it will otherwise always */ in decNumberPower()
2314 /* need to return fully-padded 1.0000 etc., but rhsint->1 */ in decNumberPower()
2315 *dac->lsu=1; /* was 0, make int 1 */ in decNumberPower()
2317 Int shift=set->digits-1; in decNumberPower()
2318 dac->digits=decShiftToMost(dac->lsu, 1, shift); in decNumberPower()
2319 dac->exponent=-shift; /* make 1.0000... */ in decNumberPower()
2328 } /* non-integer rhs */ in decNumberPower()
2332 *dac->lsu=1; /* .. */ in decNumberPower()
2336 if (decNumberIsNegative(rhs)) { /* was a **-n [hence digits>0] */ in decNumberPower()
2340 if (set->extended) { /* need to calculate 1/lhs */ in decNumberPower()
2347 if (allocinv==NULL) { /* hopeless -- abandon */ in decNumberPower()
2352 /* [inv now points to big-enough buffer or allocated storage] */ in decNumberPower()
2361 /* Raise-to-the-power loop... */ in decNumberPower()
2362 seenbit=0; /* set once a 1-bit is encountered */ in decNumberPower()
2369 /* compiler, with symptom: 5**3 -> 25, when n=n+n was used] */ in decNumberPower()
2383 /* If subset, and power was negative, reverse the kind of -erflow */ in decNumberPower()
2385 if (!set->extended && decNumberIsNegative(rhs)) { in decNumberPower()
2388 else { /* trickier -- Underflow may or may not be set */ in decNumberPower()
2394 dac->bits=(dac->bits & ~DECNEG) | bits; /* force correct sign */ in decNumberPower()
2403 if (!set->extended && /* subset math */ in decNumberPower()
2404 decNumberIsNegative(rhs)) { /* was a **-n [hence digits>0] */ in decNumberPower()
2415 if (!set->extended) decTrim(res, set, 0, &dropped); /* trailing zeros */ in decNumberPower()
2432 /* ------------------------------------------------------------------ */
2433 /* decNumberQuantize -- force exponent to requested value */
2436 /* of C (by rounding or shifting) such that the exponent (-scale) */
2445 /* C must have space for set->digits digits. */
2449 /* ------------------------------------------------------------------ */
2458 /* ------------------------------------------------------------------ */
2459 /* decNumberReduce -- remove trailing zeros */
2467 /* C must have space for set->digits digits. */
2468 /* ------------------------------------------------------------------ */
2478 decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */ in decNumberReduce()
2490 if (!set->extended) { in decNumberReduce()
2492 if (rhs->digits>set->digits) { in decNumberReduce()
2520 /* ------------------------------------------------------------------ */
2521 /* decNumberRescale -- force exponent to requested value */
2524 /* of C (by rounding or shifting) such that the exponent (-scale) */
2533 /* C must have space for set->digits digits. */
2537 /* ------------------------------------------------------------------ */
2546 /* ------------------------------------------------------------------ */
2547 /* decNumberRemainder -- divide and return remainder */
2556 /* C must have space for set->digits digits. */
2557 /* ------------------------------------------------------------------ */
2569 /* ------------------------------------------------------------------ */
2570 /* decNumberRemainderNear -- divide and return remainder from nearest */
2579 /* C must have space for set->digits digits. */
2580 /* ------------------------------------------------------------------ */
2592 /* ------------------------------------------------------------------ */
2593 /* decNumberRotate -- rotate the coefficient of a Number left/right */
2595 /* This computes C = A rot B (in base ten and rotating set->digits */
2600 /* rhs is B, the number of digits to rotate (-ve to right) */
2605 /* the exponent or the sign of A. If lhs->digits is less than */
2606 /* set->digits the coefficient is padded with zeros on the left */
2610 /* B must be an integer (q=0) and in the range -set->digits through */
2611 /* +set->digits. */
2612 /* C must have space for set->digits digits. */
2616 /* ------------------------------------------------------------------ */
2630 else if (decNumberIsInfinite(rhs) || rhs->exponent!=0) in decNumberRotate()
2636 || abs(rotate)>set->digits) /* .. or out of range */ in decNumberRotate()
2640 /* convert -ve rotate to equivalent positive rotation */ in decNumberRotate()
2641 if (rotate<0) rotate=set->digits+rotate; in decNumberRotate()
2642 if (rotate!=0 && rotate!=set->digits /* zero or full rotation */ in decNumberRotate()
2644 /* left-rotate to do; 0 < rotate < set->digits */ in decNumberRotate()
2647 Unit *msu=res->lsu+D2U(res->digits)-1; /* current msu */ in decNumberRotate()
2648 Unit *msumax=res->lsu+D2U(set->digits)-1; /* rotation msu */ in decNumberRotate()
2650 res->digits=set->digits; /* now full-length */ in decNumberRotate()
2651 msudigits=MSUDIGITS(res->digits); /* actual digits in msu */ in decNumberRotate()
2653 /* rotation here is done in-place, in three steps */ in decNumberRotate()
2654 /* 1. shift all to least up to one unit to unit-align final */ in decNumberRotate()
2683 /* Step 1: amount to shift is the partial right-rotate count */ in decNumberRotate()
2684 rotate=set->digits-rotate; /* make it right-rotate */ in decNumberRotate()
2686 shift=rotate%DECDPUN; /* left-over digits count */ in decNumberRotate()
2688 uInt save=res->lsu[0]%powers[shift]; /* save low digit(s) */ in decNumberRotate()
2689 decShiftToLeast(res->lsu, D2U(res->digits), shift); in decNumberRotate()
2690 if (shift>msudigits) { /* msumax-1 needs >0 digits */ in decNumberRotate()
2691 uInt rem=save%powers[shift-msudigits];/* split save */ in decNumberRotate()
2692 *msumax=(Unit)(save/powers[shift-msudigits]); /* and insert */ in decNumberRotate()
2693 *(msumax-1)=*(msumax-1) in decNumberRotate()
2694 +(Unit)(rem*powers[DECDPUN-(shift-msudigits)]); /* .. */ in decNumberRotate()
2697 *msumax=*msumax+(Unit)(save*powers[msudigits-shift]); /* [maybe *1] */ in decNumberRotate()
2704 /* if any, and the shift is DECDPUN-msudigits (which may be */ in decNumberRotate()
2706 shift=DECDPUN-msudigits; in decNumberRotate()
2708 uInt save=res->lsu[0]%powers[shift]; /* save low digit(s) */ in decNumberRotate()
2709 decShiftToLeast(res->lsu, units, shift); in decNumberRotate()
2715 decReverse(res->lsu+units, msumax); /* left part */ in decNumberRotate()
2716 decReverse(res->lsu, res->lsu+units-1); /* right part */ in decNumberRotate()
2717 decReverse(res->lsu, msumax); /* whole */ in decNumberRotate()
2721 res->digits=decGetDigits(res->lsu, msumax-res->lsu+1); in decNumberRotate()
2729 /* ------------------------------------------------------------------ */
2730 /* decNumberSameQuantum -- test for equal exponents */
2737 /* ------------------------------------------------------------------ */
2751 else if (lhs->exponent==rhs->exponent) ret=1; in decNumberSameQuantum()
2754 *res->lsu=ret; in decNumberSameQuantum()
2758 /* ------------------------------------------------------------------ */
2759 /* decNumberScaleB -- multiply by a power of 10 */
2769 /* C must have space for set->digits digits. */
2772 /* ------------------------------------------------------------------ */
2787 else if (decNumberIsInfinite(rhs) || rhs->exponent!=0) in decNumberScaleB()
2794 || abs(reqexp)>(2*(set->digits+set->emax))) /* .. or out of range */ in decNumberScaleB()
2799 res->exponent+=reqexp; /* adjust the exponent */ in decNumberScaleB()
2809 /* ------------------------------------------------------------------ */
2810 /* decNumberShift -- shift the coefficient of a Number left or right */
2812 /* This computes C = A << B or C = A >> -B (in base ten). */
2816 /* rhs is B, the number of digits to shift (-ve to right) */
2823 /* B must be an integer (q=0) and in the range -set->digits through */
2824 /* +set->digits. */
2825 /* C must have space for set->digits digits. */
2829 /* ------------------------------------------------------------------ */
2843 else if (decNumberIsInfinite(rhs) || rhs->exponent!=0) in decNumberShift()
2849 || abs(shift)>set->digits) /* .. or out of range */ in decNumberShift()
2855 if (shift==set->digits) { /* removing all */ in decNumberShift()
2856 *res->lsu=0; /* so place 0 */ in decNumberShift()
2857 res->digits=1; /* .. */ in decNumberShift()
2860 /* first remove leading digits if necessary */ in decNumberShift()
2861 if (res->digits+shift>set->digits) { in decNumberShift()
2862 decDecap(res, res->digits+shift-set->digits); in decNumberShift()
2863 /* that updated res->digits; may have gone to 1 (for a */ in decNumberShift()
2866 if (res->digits>1 || *res->lsu) /* if non-zero.. */ in decNumberShift()
2867 res->digits=decShiftToMost(res->lsu, res->digits, shift); in decNumberShift()
2871 if (-shift>=res->digits) { /* discarding all */ in decNumberShift()
2872 *res->lsu=0; /* so place 0 */ in decNumberShift()
2873 res->digits=1; /* .. */ in decNumberShift()
2876 decShiftToLeast(res->lsu, D2U(res->digits), -shift); in decNumberShift()
2877 res->digits-=(-shift); in decNumberShift()
2880 } /* non-0 non-Inf shift */ in decNumberShift()
2887 /* ------------------------------------------------------------------ */
2888 /* decNumberSquareRoot -- square root operator */
2896 /* C must have space for set->digits digits. */
2897 /* ------------------------------------------------------------------ */
2898 /* This uses the following varying-precision algorithm in: */
2902 /* pp229-237, ACM, September 1985. */
2904 /* The square-root is calculated using Newton's method, after which */
2934 /* p := min(2*p - 2, maxp) % p = 4,6,10, . . . , maxp */
2941 /* % of f; to ensure proper rounding, compare squares of (approx - */
2946 /* const approxsubhalf := approx - setexp(.5, -p) */
2948 /* approx := approx - setexp(.l, -p + 1) */
2950 /* const approxaddhalf := approx + setexp(.5, -p) */
2952 /* approx := approx + setexp(.l, -p + 1) */
2958 /* ------------------------------------------------------------------ */
2974 decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */ in decNumberSquareRoot()
2982 decNumber *allocbuff=NULL; /* -> allocated buff, iff allocated */ in decNumberSquareRoot()
2983 decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */ in decNumberSquareRoot()
2984 decNumber *allocbufb=NULL; /* -> allocated bufb, iff allocated */ in decNumberSquareRoot()
2990 decNumber *t=buft; /* up-to-3-digit constant or work */ in decNumberSquareRoot()
2998 if (!set->extended) { in decNumberSquareRoot()
3000 if (rhs->digits>set->digits) { in decNumberSquareRoot()
3022 /* [We would like to write: ideal=rhs->exponent>>1, but this */ in decNumberSquareRoot()
3024 ideal=(rhs->exponent&~1)/2; /* target */ in decNumberSquareRoot()
3028 decNumberCopy(res, rhs); /* could be 0 or -0 */ in decNumberSquareRoot()
3029 res->exponent=ideal; /* use the ideal [safe] */ in decNumberSquareRoot()
3030 /* use decFinish to clamp any out-of-range exponent, etc. */ in decNumberSquareRoot()
3035 /* any other -x is an oops */ in decNumberSquareRoot()
3042 /* f -- the same precision as the RHS, reduced to 0.01->0.99... */ in decNumberSquareRoot()
3043 /* a -- Hull's approximation -- precision, when assigned, is */ in decNumberSquareRoot()
3046 /* b -- intermediate temporary result (same size as a) */ in decNumberSquareRoot()
3048 workp=MAXI(set->digits+1, rhs->digits); /* actual rounding precision */ in decNumberSquareRoot()
3051 needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit); in decNumberSquareRoot()
3054 if (allocbuff==NULL) { /* hopeless -- abandon */ in decNumberSquareRoot()
3059 /* a and b both need to be able to hold a maxp-length number */ in decNumberSquareRoot()
3060 needbytes=sizeof(decNumber)+(D2U(maxp)-1)*sizeof(Unit); in decNumberSquareRoot()
3071 /* copy rhs -> f, save exponent, and reduce so 0.1 <= f < 1 */ in decNumberSquareRoot()
3073 exp=f->exponent+f->digits; /* adjusted to Hull rules */ in decNumberSquareRoot()
3074 f->exponent=-(f->digits); /* to range */ in decNumberSquareRoot()
3084 t->bits=0; t->digits=3; in decNumberSquareRoot()
3085 a->bits=0; a->digits=3; in decNumberSquareRoot()
3088 t->exponent=-3; in decNumberSquareRoot()
3089 a->exponent=-3; in decNumberSquareRoot()
3091 t->lsu[0]=259; in decNumberSquareRoot()
3092 a->lsu[0]=819; in decNumberSquareRoot()
3094 t->lsu[0]=59; t->lsu[1]=2; in decNumberSquareRoot()
3095 a->lsu[0]=19; a->lsu[1]=8; in decNumberSquareRoot()
3097 t->lsu[0]=9; t->lsu[1]=5; t->lsu[2]=2; in decNumberSquareRoot()
3098 a->lsu[0]=9; a->lsu[1]=1; a->lsu[2]=8; in decNumberSquareRoot()
3103 f->exponent--; /* f=f/10 */ in decNumberSquareRoot()
3105 t->exponent=-4; in decNumberSquareRoot()
3106 a->exponent=-2; in decNumberSquareRoot()
3108 t->lsu[0]=819; in decNumberSquareRoot()
3109 a->lsu[0]=259; in decNumberSquareRoot()
3111 t->lsu[0]=19; t->lsu[1]=8; in decNumberSquareRoot()
3112 a->lsu[0]=59; a->lsu[1]=2; in decNumberSquareRoot()
3114 t->lsu[0]=9; t->lsu[1]=1; t->lsu[2]=8; in decNumberSquareRoot()
3115 a->lsu[0]=9; a->lsu[1]=5; a->lsu[2]=2; in decNumberSquareRoot()
3126 t->lsu[0]=5; /* .. */ in decNumberSquareRoot()
3127 t->exponent=-1; /* .. */ in decNumberSquareRoot()
3130 /* set p to min(2*p - 2, maxp) [hence 3; or: 4, 6, 10, ... , maxp] */ in decNumberSquareRoot()
3131 workset.digits=workset.digits*2-2; in decNumberSquareRoot()
3138 if (a->digits==maxp) break; /* have required digits */ in decNumberSquareRoot()
3147 a->exponent+=exp/2; /* set correct exponent */ in decNumberSquareRoot()
3154 /* Overflow was possible if the input exponent was out-of-range, */ in decNumberSquareRoot()
3157 status=rstatus; /* use the status as-is */ in decNumberSquareRoot()
3166 a->exponent-=exp/2; /* back to 0.1->1 */ in decNumberSquareRoot()
3170 /* squares of (a - l/2 ulp) and (a + l/2 ulp) with f. */ in decNumberSquareRoot()
3171 /* Here workset.digits=maxp and t=0.5, and a->digits determines */ in decNumberSquareRoot()
3173 workset.digits--; /* maxp-1 is OK now */ in decNumberSquareRoot()
3174 t->exponent=-a->digits-1; /* make 0.5 ulp */ in decNumberSquareRoot()
3175 decAddOp(b, a, t, &workset, DECNEG, &ignore); /* b = a - 0.5 ulp */ in decNumberSquareRoot()
3181 t->exponent++; /* make 1.0 ulp */ in decNumberSquareRoot()
3182 t->lsu[0]=1; /* .. */ in decNumberSquareRoot()
3183 decAddOp(a, a, t, &workset, DECNEG, &ignore); /* a = a - 1 ulp */ in decNumberSquareRoot()
3185 approxset.emin-=exp/2; /* adjust to match a */ in decNumberSquareRoot()
3186 approxset.emax-=exp/2; in decNumberSquareRoot()
3195 t->exponent++; /* make 1.0 ulp */ in decNumberSquareRoot()
3196 t->lsu[0]=1; /* .. */ in decNumberSquareRoot()
3199 approxset.emin-=exp/2; /* adjust to match a */ in decNumberSquareRoot()
3200 approxset.emax-=exp/2; in decNumberSquareRoot()
3208 a->exponent+=exp/2; /* set correct exponent */ in decNumberSquareRoot()
3219 if (b->digits*2-1 > workp && !set->clamp) { /* cannot fit */ in decNumberSquareRoot()
3234 Int todrop=ideal-a->exponent; /* most that can be dropped */ in decNumberSquareRoot()
3242 decShiftToLeast(a->lsu, D2U(a->digits), todrop); in decNumberSquareRoot()
3243 a->exponent+=todrop; /* maintain numerical value */ in decNumberSquareRoot()
3244 a->digits-=todrop; /* new length */ in decNumberSquareRoot()
3251 /* double-check Underflow, as perhaps the result could not have */ in decNumberSquareRoot()
3254 Int ae=rhs->exponent+rhs->digits-1; /* adjusted exponent */ in decNumberSquareRoot()
3257 if (ae>=set->emin*2) status&=~(DEC_Subnormal|DEC_Underflow); in decNumberSquareRoot()
3259 if (ae>=set->emin*2) status&=~DEC_Underflow; in decNumberSquareRoot()
3281 /* ------------------------------------------------------------------ */
3282 /* decNumberSubtract -- subtract two Numbers */
3284 /* This computes C = A - B */
3286 /* res is C, the result. C may be A and/or B (e.g., X=X-X) */
3291 /* C must have space for set->digits digits. */
3292 /* ------------------------------------------------------------------ */
3305 /* ------------------------------------------------------------------ */
3306 /* decNumberToIntegralExact -- round-to-integral-value with InExact */
3307 /* decNumberToIntegralValue -- round-to-integral-value */
3317 /* rescale(rhs, 0) if rhs->exponent is <0. */
3325 /* ------------------------------------------------------------------ */
3343 if (rhs->exponent>=0) return decNumberCopy(res, rhs); in decNumberToIntegralExact()
3346 workset.digits=rhs->digits; /* no length rounding */ in decNumberToIntegralExact()
3363 set->status|=workset.status&DEC_Invalid_operation; in decNumberToIntegralValue()
3367 /* ------------------------------------------------------------------ */
3368 /* decNumberXor -- XOR two Numbers, digitwise */
3377 /* C must have space for set->digits digits. */
3381 /* ------------------------------------------------------------------ */
3384 const Unit *ua, *ub; /* -> operands */ in decNumberXor()
3385 const Unit *msua, *msub; /* -> operand msus */ in decNumberXor()
3386 Unit *uc, *msuc; /* -> result and its msu */ in decNumberXor()
3392 if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs) in decNumberXor()
3393 || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) { in decNumberXor()
3398 ua=lhs->lsu; /* bottom-up */ in decNumberXor()
3399 ub=rhs->lsu; /* .. */ in decNumberXor()
3400 uc=res->lsu; /* .. */ in decNumberXor()
3401 msua=ua+D2U(lhs->digits)-1; /* -> msu of lhs */ in decNumberXor()
3402 msub=ub+D2U(rhs->digits)-1; /* -> msu of rhs */ in decNumberXor()
3403 msuc=uc+D2U(set->digits)-1; /* -> msu of result */ in decNumberXor()
3404 msudigs=MSUDIGITS(set->digits); /* [faster than remainder] */ in decNumberXor()
3425 if (uc==msuc && i==msudigs-1) break; /* just did final digit */ in decNumberXor()
3427 } /* non-zero */ in decNumberXor()
3429 /* [here uc-1 is the msu of the result] */ in decNumberXor()
3430 res->digits=decGetDigits(res->lsu, uc-res->lsu); in decNumberXor()
3431 res->exponent=0; /* integer */ in decNumberXor()
3432 res->bits=0; /* sign=0 */ in decNumberXor()
3441 /* ------------------------------------------------------------------ */
3442 /* decNumberClass -- return the decClass of a decNumber */
3443 /* dn -- the decNumber to test */
3444 /* set -- the context to use for Emin */
3446 /* ------------------------------------------------------------------ */
3469 /* ------------------------------------------------------------------ */
3470 /* decNumberClassToString -- convert decClass to a string */
3474 /* ------------------------------------------------------------------ */
3489 /* ------------------------------------------------------------------ */
3490 /* decNumberCopy -- copy a number */
3496 /* (dest==src is allowed and is a no-op) */
3499 /* ------------------------------------------------------------------ */
3512 dest->bits=src->bits; in decNumberCopy()
3513 dest->exponent=src->exponent; in decNumberCopy()
3514 dest->digits=src->digits; in decNumberCopy()
3515 dest->lsu[0]=src->lsu[0]; in decNumberCopy()
3516 if (src->digits>DECDPUN) { /* more Units to come */ in decNumberCopy()
3521 d=dest->lsu+1; /* -> first destination */ in decNumberCopy()
3522 smsup=src->lsu+D2U(src->digits); /* -> source msu+1 */ in decNumberCopy()
3523 for (s=src->lsu+1; s<smsup; s++, d++) *d=*s; in decNumberCopy()
3528 /* ------------------------------------------------------------------ */
3529 /* decNumberCopyAbs -- quiet absolute value operator */
3536 /* C must have space for set->digits digits. */
3539 /* ------------------------------------------------------------------ */
3545 res->bits&=~DECNEG; /* turn off sign */ in decNumberCopyAbs()
3549 /* ------------------------------------------------------------------ */
3550 /* decNumberCopyNegate -- quiet negate value operator */
3557 /* C must have space for set->digits digits. */
3560 /* ------------------------------------------------------------------ */
3566 res->bits^=DECNEG; /* invert the sign */ in decNumberCopyNegate()
3570 /* ------------------------------------------------------------------ */
3571 /* decNumberCopySign -- quiet copy and set sign operator */
3579 /* C must have space for set->digits digits. */
3581 /* ------------------------------------------------------------------ */
3588 sign=rhs->bits & DECNEG; /* save sign bit */ in decNumberCopySign()
3590 res->bits&=~DECNEG; /* clear the sign */ in decNumberCopySign()
3591 res->bits|=sign; /* set from rhs */ in decNumberCopySign()
3595 /* ------------------------------------------------------------------ */
3596 /* decNumberGetBCD -- get the coefficient in BCD8 */
3598 /* bcd is the uInt array that will receive dn->digits BCD bytes, */
3599 /* most-significant at offset 0 */
3602 /* bcd must have at least dn->digits bytes. No error is possible; if */
3604 /* ------------------------------------------------------------------ */
3606 uByte *ub=bcd+dn->digits-1; /* -> lsd */ in decNumberGetBCD()
3607 const Unit *up=dn->lsu; /* Unit pointer, -> lsu */ in decNumberGetBCD()
3610 for (; ub>=bcd; ub--, up++) *ub=*up; in decNumberGetBCD()
3614 for (; ub>=bcd; ub--) { in decNumberGetBCD()
3617 cut--; in decNumberGetBCD()
3627 /* ------------------------------------------------------------------ */
3628 /* decNumberSetBCD -- set (replace) the coefficient from BCD8 */
3630 /* bcd is the uInt array that will source n BCD bytes, most- */
3638 /* ------------------------------------------------------------------ */
3640 Unit *up = dn->lsu + D2U(n) - 1; /* -> msu [target pointer] */ in decNumberSetBCD()
3641 const uByte *ub=bcd; /* -> source msd */ in decNumberSetBCD()
3644 for (; ub<bcd+n; ub++, up--) *up=*ub; in decNumberSetBCD()
3648 for (;up>=dn->lsu; up--) { /* each Unit from msu */ in decNumberSetBCD()
3650 for (; cut>0; ub++, cut--) *up=X10(*up)+*ub; in decNumberSetBCD()
3654 dn->digits=n; /* set digit count */ in decNumberSetBCD()
3658 /* ------------------------------------------------------------------ */
3659 /* decNumberIsNormal -- test normality of a decNumber */
3663 /* ------------------------------------------------------------------ */
3671 if (decNumberIsZero(dn)) return 0; /* not non-zero */ in decNumberIsNormal()
3673 ae=dn->exponent+dn->digits-1; /* adjusted exponent */ in decNumberIsNormal()
3674 if (ae<set->emin) return 0; /* is subnormal */ in decNumberIsNormal()
3678 /* ------------------------------------------------------------------ */
3679 /* decNumberIsSubnormal -- test subnormality of a decNumber */
3682 /* returns 1 if |dn| is finite, non-zero, and <Nmin, 0 otherwise */
3683 /* ------------------------------------------------------------------ */
3691 if (decNumberIsZero(dn)) return 0; /* not non-zero */ in decNumberIsSubnormal()
3693 ae=dn->exponent+dn->digits-1; /* adjusted exponent */ in decNumberIsSubnormal()
3694 if (ae<set->emin) return 1; /* is subnormal */ in decNumberIsSubnormal()
3698 /* ------------------------------------------------------------------ */
3699 /* decNumberTrim -- remove insignificant zeros */
3706 /* ------------------------------------------------------------------ */
3717 /* ------------------------------------------------------------------ */
3718 /* decNumberVersion -- return the name and version of this module */
3721 /* ------------------------------------------------------------------ */
3726 /* ------------------------------------------------------------------ */
3727 /* decNumberZero -- set a number to 0 */
3733 /* ------------------------------------------------------------------ */
3741 dn->bits=0; in decNumberZero()
3742 dn->exponent=0; in decNumberZero()
3743 dn->digits=1; in decNumberZero()
3744 dn->lsu[0]=0; in decNumberZero()
3752 /* ------------------------------------------------------------------ */
3753 /* decToString -- lay out a number into a string */
3759 /* string must be at least dn->digits+14 characters long */
3762 /* Note that this routine can generate a -0 or 0.000. These are */
3763 /* never generated in subset to-number or arithmetic, but can occur */
3764 /* in non-subset arithmetic (e.g., -1*0 or 1.234-1.234). */
3765 /* ------------------------------------------------------------------ */
3769 Int exp=dn->exponent; /* local copy */ in decToString()
3770 Int e; /* E-part value */ in decToString()
3774 const Unit *up=dn->lsu+D2U(dn->digits)-1; /* -> msu [input pointer] */ in decToString()
3784 *c='-'; in decToString()
3787 if (dn->bits&DECSPECIAL) { /* Is a special value */ in decToString()
3793 if (dn->bits&DECSNAN) { /* signalling NaN */ in decToString()
3799 /* if not a clean non-zero coefficient, that's all there is in a */ in decToString()
3801 if (exp!=0 || (*dn->lsu==0 && dn->digits==1)) return; in decToString()
3806 cut=MSUDIGITS(dn->digits); /* [faster than remainder] */ in decToString()
3807 cut--; /* power of ten for digit */ in decToString()
3810 for (;up>=dn->lsu; up--) { /* each Unit from msu */ in decToString()
3812 for (; cut>=0; c++, cut--) TODIGIT(u, cut, c, pow); in decToString()
3813 cut=DECDPUN-1; /* next Unit has all digits */ in decToString()
3818 /* non-0 exponent -- assume plain form */ in decToString()
3819 pre=dn->digits+exp; /* digits before '.' */ in decToString()
3821 if ((exp>0) || (pre<-5)) { /* need exponential form */ in decToString()
3822 e=exp+dn->digits-1; /* calculate E value */ in decToString()
3829 adj=(-e)%3; in decToString()
3830 if (adj!=0) adj=3-adj; in decToString()
3835 e=e-adj; in decToString()
3843 pre=-(2-adj); in decToString()
3853 for (; pre>0; pre--, c++, cut--) { in decToString()
3855 if (up==dn->lsu) break; /* out of input digits (pre>digits) */ in decToString()
3856 up--; in decToString()
3857 cut=DECDPUN-1; in decToString()
3862 if (n<dn->digits) { /* more to come, after '.' */ in decToString()
3864 for (;; c++, cut--) { in decToString()
3866 if (up==dn->lsu) break; /* out of input digits */ in decToString()
3867 up--; in decToString()
3868 cut=DECDPUN-1; in decToString()
3874 else for (; pre>0; pre--, c++) *c='0'; /* 0 padding (for engineering) needed */ in decToString()
3880 for (; ; c++, cut--) { in decToString()
3882 if (up==dn->lsu) break; /* out of input digits */ in decToString()
3883 up--; in decToString()
3884 cut=DECDPUN-1; in decToString()
3891 /* Finally add the E-part, if needed. It will never be 0, has a in decToString()
3892 base maximum and minimum of +999999999 through -999999999, but in decToString()
3893 could range down to -1999999998 for anormal numbers */ in decToString()
3895 Flag had=0; /* 1=had non-zero */ in decToString()
3900 *(c-1)='-'; /* oops, need - */ in decToString()
3901 u=-e; /* uInt, please */ in decToString()
3904 for (cut=9; cut>=0; cut--) { in decToString()
3907 had=1; /* had non-0 */ in decToString()
3915 /* ------------------------------------------------------------------ */
3916 /* decAddOp -- add/subtract operation */
3927 /* C must have space for set->digits digits. */
3929 /* ------------------------------------------------------------------ */
3932 /* -- a digits+1 calculation is needed because the numbers are */
3933 /* unaligned and span more than set->digits digits */
3934 /* -- a carry to digits+1 digits looks possible */
3935 /* -- C is the same as A or B, and the result would destructively */
3946 /* Addition, especially x=x+1, is speed-critical. */
3948 /* calls from higher-level functions (notably exp). */
3949 /* ------------------------------------------------------------------ */
3954 decNumber *alloclhs=NULL; /* non-NULL if rounded lhs allocated */ in decAddOp()
3962 Flag diffsign; /* non-0 if arguments have different sign */ in decAddOp()
3967 Unit *allocacc=NULL; /* -> allocated acc buffer, iff allocated */ in decAddOp()
3968 Int reqdigits=set->digits; /* local copy; requested DIGITS */ in decAddOp()
3977 if (!set->extended) { in decAddOp()
3979 if (lhs->digits>reqdigits) { in decAddOp()
3984 if (rhs->digits>reqdigits) { in decAddOp()
3994 diffsign=(Flag)((lhs->bits^rhs->bits^negate)&DECNEG); in decAddOp()
4007 bits=lhs->bits & DECNEG; /* get sign from LHS */ in decAddOp()
4009 else bits=(rhs->bits^negate) & DECNEG;/* RHS must be Infinity */ in decAddOp()
4012 res->bits=bits; /* set +/- infinity */ in decAddOp()
4017 /* Quick exit for add 0s; return the non-0, modified as need be */ in decAddOp()
4020 Int lexp=lhs->exponent; /* save in case LHS==RES */ in decAddOp()
4021 bits=lhs->bits; /* .. */ in decAddOp()
4024 res->bits^=negate; /* flip if rhs was negated */ in decAddOp()
4026 if (set->extended) { /* exponents on zeros count */ in decAddOp()
4029 adjust=lexp-res->exponent; /* adjustment needed [if -ve] */ in decAddOp()
4031 if (adjust<0) res->exponent=lexp; /* set exponent */ in decAddOp()
4032 /* 0-0 gives +0 unless rounding to -infinity, and -0-0 gives -0 */ in decAddOp()
4034 if (set->round!=DEC_ROUND_FLOOR) res->bits=0; in decAddOp()
4035 else res->bits=DECNEG; /* preserve 0 sign */ in decAddOp()
4038 else { /* non-0 res */ in decAddOp()
4039 if (adjust<0) { /* 0-padding needed */ in decAddOp()
4040 if ((res->digits-adjust)>set->digits) { in decAddOp()
4041 adjust=res->digits-set->digits; /* to fit exactly */ in decAddOp()
4044 res->digits=decShiftToMost(res->lsu, res->digits, -adjust); in decAddOp()
4045 res->exponent+=adjust; /* set the exponent. */ in decAddOp()
4047 } /* non-0 res */ in decAddOp()
4054 if (ISZERO(rhs)) { /* [lhs is non-zero] */ in decAddOp()
4056 Int rexp=rhs->exponent; /* save in case RHS==RES */ in decAddOp()
4057 bits=rhs->bits; /* be clean */ in decAddOp()
4061 if (set->extended) { /* exponents on zeros count */ in decAddOp()
4064 /* [0-0 case handled above] */ in decAddOp()
4065 adjust=rexp-res->exponent; /* adjustment needed [if -ve] */ in decAddOp()
4066 if (adjust<0) { /* 0-padding needed */ in decAddOp()
4067 if ((res->digits-adjust)>set->digits) { in decAddOp()
4068 adjust=res->digits-set->digits; /* to fit exactly */ in decAddOp()
4071 res->digits=decShiftToMost(res->lsu, res->digits, -adjust); in decAddOp()
4072 res->exponent+=adjust; /* set the exponent. */ in decAddOp()
4081 /* (notably 0-0) have already been handled] */ in decAddOp()
4084 padding=rhs->exponent-lhs->exponent; in decAddOp()
4090 && rhs->digits<=DECDPUN in decAddOp()
4091 && rhs->exponent>=set->emin /* [some normals drop through] */ in decAddOp()
4092 && rhs->exponent<=set->emax-set->digits+1 /* [could clamp] */ in decAddOp()
4093 && rhs->digits<=reqdigits in decAddOp()
4094 && lhs->digits<=reqdigits) { in decAddOp()
4095 Int partial=*lhs->lsu; in decAddOp()
4097 partial+=*rhs->lsu; in decAddOp()
4099 && (lhs->digits>=DECDPUN || /* .. and no digits-count change */ in decAddOp()
4100 partial<(Int)powers[lhs->digits])) { /* .. */ in decAddOp()
4102 *res->lsu=(Unit)partial; /* [copy could have overwritten RHS] */ in decAddOp()
4108 partial-=*rhs->lsu; in decAddOp()
4109 if (partial>0) { /* no borrow needed, and non-0 result */ in decAddOp()
4111 *res->lsu=(Unit)partial; in decAddOp()
4113 res->digits=decGetDigits(res->lsu, D2U(res->digits)); in decAddOp()
4124 /* other) padding with up to DIGITS-1 trailing zeros may be */ in decAddOp()
4128 bits=lhs->bits; /* assume sign is that of LHS */ in decAddOp()
4139 padding=-padding; /* will be +ve */ in decAddOp()
4140 bits=(uByte)(rhs->bits^negate); /* assumed sign is now that of RHS */ in decAddOp()
4148 if (rhs->digits+padding > lhs->digits+reqdigits+1) { in decAddOp()
4151 Int shift=reqdigits-rhs->digits; /* left shift needed */ in decAddOp()
4153 if (diffsign) residue=-residue; /* signs differ */ in decAddOp()
4158 res->digits=decShiftToMost(res->lsu, res->digits, shift); in decAddOp()
4159 res->exponent-=shift; /* adjust the exponent. */ in decAddOp()
4162 if (!swapped) res->bits^=negate; in decAddOp()
4167 rhsshift=D2U(padding+1)-1; /* this much by Unit shift .. */ in decAddOp()
4168 mult=powers[padding-(rhsshift*DECDPUN)]; /* .. this by multiplication */ in decAddOp()
4171 if (diffsign) mult=-mult; /* signs differ */ in decAddOp()
4174 maxdigits=rhs->digits+padding; /* virtual length of RHS */ in decAddOp()
4175 if (lhs->digits>maxdigits) maxdigits=lhs->digits; in decAddOp()
4179 acc=res->lsu; /* assume add direct to result */ in decAddOp()
4192 if (allocacc==NULL) { /* hopeless -- abandon */ in decAddOp()
4199 res->bits=(uByte)(bits&DECNEG); /* it's now safe to overwrite.. */ in decAddOp()
4200 res->exponent=lhs->exponent; /* .. operands (even if aliased) */ in decAddOp()
4203 decDumpAr('A', lhs->lsu, D2U(lhs->digits)); in decAddOp()
4204 decDumpAr('B', rhs->lsu, D2U(rhs->digits)); in decAddOp()
4208 /* add [A+B*m] or subtract [A+B*(-m)] */ in decAddOp()
4209 res->digits=decUnitAddSub(lhs->lsu, D2U(lhs->digits), in decAddOp()
4210 rhs->lsu, D2U(rhs->digits), in decAddOp()
4212 *DECDPUN; /* [units -> digits] */ in decAddOp()
4213 if (res->digits<0) { /* borrowed... */ in decAddOp()
4214 res->digits=-res->digits; in decAddOp()
4215 res->bits^=DECNEG; /* flip the sign */ in decAddOp()
4218 decDumpAr('+', acc, D2U(res->digits)); in decAddOp()
4225 if (acc!=res->lsu) { in decAddOp()
4227 if (set->extended) { /* round from first significant digit */ in decAddOp()
4229 /* remove leading zeros that were added due to rounding up to */ in decAddOp()
4230 /* integral Units -- before the test for rounding. */ in decAddOp()
4231 if (res->digits>reqdigits) in decAddOp()
4232 res->digits=decGetDigits(acc, D2U(res->digits)); in decAddOp()
4233 decSetCoeff(res, set, acc, res->digits, &residue, status); in decAddOp()
4239 /* negative multiple (-10, -100...) and the top digit(s) become */ in decAddOp()
4242 if (res->digits<maxdigits) { in decAddOp()
4243 *(acc+D2U(res->digits))=0; /* ensure leading 0 is there */ in decAddOp()
4244 res->digits=maxdigits; in decAddOp()
4247 /* remove leading zeros that added due to rounding up to */ in decAddOp()
4250 if (res->digits>reqdigits) { in decAddOp()
4251 res->digits=decGetDigits(acc, D2U(res->digits)); in decAddOp()
4252 if (res->digits<maxdigits) res->digits=maxdigits; in decAddOp()
4255 decSetCoeff(res, set, acc, res->digits, &residue, status); in decAddOp()
4267 res->digits=decGetDigits(res->lsu, D2U(res->digits)); in decAddOp()
4274 /* except round toward -Infinity, in which mode that sign shall be */ in decAddOp()
4275 /* '-'." [Subset zeros also never have '-', set by decFinish.] */ in decAddOp()
4278 && set->extended in decAddOp()
4281 if (set->round==DEC_ROUND_FLOOR) res->bits|=DECNEG; /* sign - */ in decAddOp()
4282 else res->bits&=~DECNEG; /* sign + */ in decAddOp()
4294 /* ------------------------------------------------------------------ */
4295 /* decDivideOp -- division operation */
4309 /* C must have space for set->digits digits. */
4311 /* ------------------------------------------------------------------ */
4313 /* 1981 S/370 implementation, that is, non-restoring long division */
4314 /* with bi-unit (rather than bi-digit) estimation for each unit */
4321 /* Exp =Exp1 - Exp2 */
4322 /* Exp =Exp +len(var1) -len(var2) */
4324 /* Pad accumulator (Var1) to double-length with 0's (pad1) */
4338 /* If same then tops2=msu2pair -- {units 1&2 of var2} */
4339 /* else tops2=msu2plus -- {0, unit 1 of var2} */
4341 /* mult=tops1/tops2 -- Good and safe guess at divisor */
4350 /* exp=exp-1 */
4352 /* exp=exp+1 -- set the proper exponent */
4356 /* ------------------------------------------------------------------ */
4362 /* for calls from higher-level functions (notably exp). */
4363 /* ------------------------------------------------------------------ */
4368 decNumber *alloclhs=NULL; /* non-NULL if rounded lhs allocated */ in decDivideOp()
4372 Unit *acc=accbuff; /* -> accumulator array for result */ in decDivideOp()
4373 Unit *allocacc=NULL; /* -> allocated buffer, iff allocated */ in decDivideOp()
4374 Unit *accnext; /* -> where next digit will go */ in decDivideOp()
4380 Unit *var1=varbuff; /* -> var1 array for long subtraction */ in decDivideOp()
4381 Unit *varalloc=NULL; /* -> allocated buffer, iff used */ in decDivideOp()
4382 Unit *msu1; /* -> msu of var1 */ in decDivideOp()
4384 const Unit *var2; /* -> var2 array */ in decDivideOp()
4385 const Unit *msu2; /* -> msu of var2 */ in decDivideOp()
4396 Int reqdigits=set->digits; /* requested DIGITS */ in decDivideOp()
4414 if (!set->extended) { in decDivideOp()
4416 if (lhs->digits>reqdigits) { in decDivideOp()
4421 if (rhs->digits>reqdigits) { in decDivideOp()
4430 bits=(lhs->bits^rhs->bits)&DECNEG; /* assumed sign for divisions */ in decDivideOp()
4447 res->bits=bits|DECINF; /* set +/- infinity */ in decDivideOp()
4458 res->bits=bits; /* set +/- zero */ in decDivideOp()
4462 res->exponent=set->emin-set->digits+1; in decDivideOp()
4482 res->bits=bits|DECINF; /* .. is +/- Infinity */ in decDivideOp()
4490 if (!set->extended) decNumberZero(res); in decDivideOp()
4495 exponent=lhs->exponent-rhs->exponent; /* ideal exponent */ in decDivideOp()
4497 res->bits=bits; /* sign as computed */ in decDivideOp()
4498 res->exponent=exponent; /* exponent, too */ in decDivideOp()
4503 res->bits=bits; /* sign as computed */ in decDivideOp()
4506 exponent=rhs->exponent; /* [save in case overwrite] */ in decDivideOp()
4508 if (exponent<res->exponent) res->exponent=exponent; /* use lower */ in decDivideOp()
4519 exponent=(lhs->exponent+lhs->digits)-(rhs->exponent+rhs->digits); in decDivideOp()
4521 /* If the working exponent is -ve, then some quick exits are */ in decDivideOp()
4523 /* [for REMNEAR, it needs to be < -1, as -0.5 could need work] */ in decDivideOp()
4528 if (set->extended) in decDivideOp()
4530 res->bits=bits; /* set +/- zero */ in decDivideOp()
4534 if (lhs->exponent<=rhs->exponent) { in decDivideOp()
4535 if (op&REMAINDER || exponent<-1) { in decDivideOp()
4537 /* clone of] lhs (r = x - 0*y) */ in decDivideOp()
4555 if (allocacc==NULL) { /* hopeless -- abandon */ in decDivideOp()
4565 /* (rhs->digits+reqdigits-1) -- to allow full slide to right */ in decDivideOp()
4566 /* or (lhs->digits) -- to allow for long lhs */ in decDivideOp()
4568 /* +1 -- for rounding of slide to right */ in decDivideOp()
4569 /* +1 -- for leading 0s */ in decDivideOp()
4570 /* +1 -- for pre-adjust if a remainder or DIVIDEINT */ in decDivideOp()
4572 maxdigits=rhs->digits+reqdigits-1; in decDivideOp()
4573 if (lhs->digits>maxdigits) maxdigits=lhs->digits; in decDivideOp()
4580 if (varalloc==NULL) { /* hopeless -- abandon */ in decDivideOp()
4592 msu1=var1+var1units-1; /* msu of var1 */ in decDivideOp()
4593 source=lhs->lsu+D2U(lhs->digits)-1; /* msu of input array */ in decDivideOp()
4594 for (target=msu1; source>=lhs->lsu; source--, target--) *target=*source; in decDivideOp()
4595 for (; target>=var1; target--) *target=0; in decDivideOp()
4597 /* rhs (var2) is left-aligned with var1 at the start */ in decDivideOp()
4599 var2units=D2U(rhs->digits); /* rhs actual length (units) */ in decDivideOp()
4600 var2=rhs->lsu; /* -> rhs array */ in decDivideOp()
4601 msu2=var2+var2units-1; /* -> msu of var2 [never changes] */ in decDivideOp()
4609 msu2pair+=*(msu2-1); /* .. */ in decDivideOp()
4615 /* both left-aligned. Adjust the exponent to compensate: add the */ in decDivideOp()
4618 /* lead1=DECDPUN-digits1, and similarly for lead2.] */ in decDivideOp()
4619 for (pow=&powers[1]; *msu1>=*pow; pow++) exponent--; in decDivideOp()
4623 /* the result will be Unit-aligned. To do this, shift the var1 */ in decDivideOp()
4630 var1initpad=(var1units-D2U(lhs->digits))*DECDPUN; in decDivideOp()
4632 if (exponent<0) cut=-exponent; in decDivideOp()
4633 else cut=DECDPUN-exponent%DECDPUN; in decDivideOp()
4636 var1initpad-=cut; /* .. and reduce padding */ in decDivideOp()
4637 /* clean any most-significant units which were just emptied */ in decDivideOp()
4638 for (u=msu1; cut>=DECDPUN; cut-=DECDPUN, u--) *u=0; in decDivideOp()
4641 maxexponent=lhs->exponent-rhs->exponent; /* save */ in decDivideOp()
4645 var2ulen--; /* shift down */ in decDivideOp()
4646 exponent-=DECDPUN; /* update the exponent */ in decDivideOp()
4650 /* ---- start the long-division loops ------------------------------ */ in decDivideOp()
4653 accnext=acc+acclength-1; /* -> msu of acc [NB: allows digits+1] */ in decDivideOp()
4658 /* strip leading zero units [from either pre-adjust or from */ in decDivideOp()
4660 for (; *msu1==0 && msu1>var1; msu1--) var1units--; in decDivideOp()
4663 if (var1units==var2ulen) { /* unit-by-unit compare needed */ in decDivideOp()
4667 pv2=msu2; /* -> msu */ in decDivideOp()
4668 for (pv1=msu1; ; pv1--, pv2--) { in decDivideOp()
4669 /* v1=*pv1 -- always OK */ in decDivideOp()
4687 /* Estimate the multiplier (there's always a msu1-1)... */ in decDivideOp()
4689 mult=(Int)(((eInt)*msu1*(DECDPUNMAX+1)+*(msu1-1))/msu2pair); in decDivideOp()
4694 mult=(Int)(((eInt)*msu1*(DECDPUNMAX+1)+*(msu1-1))/msu2plus); in decDivideOp()
4699 /* subtract var1-var2, into var1; only the overlap needs */ in decDivideOp()
4700 /* processing, as this is an in-place calculation */ in decDivideOp()
4701 shift=var2ulen-var2units; in decDivideOp()
4703 decDumpAr('1', &var1[shift], var1units-shift); in decDivideOp()
4705 printf("m=%ld\n", -mult); in decDivideOp()
4707 decUnitAddSub(&var1[shift], var1units-shift, in decDivideOp()
4709 &var1[shift], -mult); in decDivideOp()
4711 decDumpAr('#', &var1[shift], var1units-shift); in decDivideOp()
4719 if (accunits!=0 || thisunit!=0) { /* is first or non-zero */ in decDivideOp()
4728 accnext--; /* ready for next */ in decDivideOp()
4743 /* to get here, var1 is less than var2, so divide var2 by the per- */ in decDivideOp()
4745 var2ulen--; /* shift down */ in decDivideOp()
4746 exponent-=DECDPUN; /* update the exponent */ in decDivideOp()
4749 /* ---- division is complete --------------------------------------- */ in decDivideOp()
4760 /* accnext now -> lowest unit of result */ in decDivideOp()
4768 /* There will be at most DECDPUN-1, from the final multiply, */ in decDivideOp()
4769 /* and then only if the result is non-0 (and even) and the */ in decDivideOp()
4779 if ((lsu-QUOT10(lsu, drop+1) in decDivideOp()
4780 *powers[drop+1])!=0) break; /* found non-0 digit */ in decDivideOp()
4782 if (lsu%powers[drop+1]!=0) break; /* found non-0 digit */ in decDivideOp()
4810 bits=lhs->bits; /* remainder sign is always as lhs */ in decDivideOp()
4815 Int exp=lhs->exponent; /* save min(exponents) */ in decDivideOp()
4816 if (rhs->exponent<exp) exp=rhs->exponent; in decDivideOp()
4819 if (set->extended) in decDivideOp()
4821 res->exponent=exp; /* .. with proper exponent */ in decDivideOp()
4822 res->bits=(uByte)(bits&DECNEG); /* [cleaned] */ in decDivideOp()
4835 postshift=var1initpad+exponent-lhs->exponent+rhs->exponent; in decDivideOp()
4845 exponent=lhs->exponent; /* exponent is smaller of lhs & rhs */ in decDivideOp()
4846 if (rhs->exponent<exponent) exponent=rhs->exponent; in decDivideOp()
4850 /* the integer was odd then the result should be rem-rhs. */ in decDivideOp()
4864 compare=decUnitCompare(accnext, tarunits, rhs->lsu, D2U(rhs->digits), in decDivideOp()
4865 rhs->exponent-exponent); in decDivideOp()
4877 *(up-1)+=DIV_ROUND_UP(DECDPUNMAX, 2); in decDivideOp()
4883 /* This is effectively causing round-up of the quotient, */ in decDivideOp()
4885 /* nines, it would overflow and hence division-impossible */ in decDivideOp()
4891 if (*up!=DECDPUNMAX) break;/* non-nines */ in decDivideOp()
4894 if (*up==powers[quotdigits]-1) allnines=1; in decDivideOp()
4897 quotdigits-=DECDPUN; /* checked those digits */ in decDivideOp()
4904 /* rem-rhs is needed; the sign will invert. Again, var1 */ in decDivideOp()
4906 exp=rhs->exponent-exponent; /* RHS padding needed */ in decDivideOp()
4910 /* subtract [A+B*(-m)]; the result will always be negative */ in decDivideOp()
4911 accunits=-decUnitAddSub(accnext, accunits, in decDivideOp()
4912 rhs->lsu, D2U(rhs->digits), in decDivideOp()
4913 expunits, accnext, -(Int)powers[exprem]); in decDivideOp()
4924 res->exponent=exponent; in decDivideOp()
4925 res->bits=(uByte)(bits&DECNEG); /* [cleaned] */ in decDivideOp()
4934 if (!set->extended && (op==DIVIDE)) decTrim(res, set, 0, &dropped); in decDivideOp()
4947 /* ------------------------------------------------------------------ */
4948 /* decMultiplyOp -- multiplication operation */
4958 /* C must have space for set->digits digits. */
4960 /* ------------------------------------------------------------------ */
4965 /* There are two major paths here: the general-purpose ('old code') */
4967 /* which is used if 64-bit ints are available, DECDPUN<=4, and more */
4970 /* The fastpath version lumps units together into 8-digit or 9-digit */
4972 /* 64-bit divisions. The chunks are then broken apart again into */
4974 /* fastpath can speed up some 16-digit operations by 10x (and much */
4975 /* more for higher-precision calculations). */
4982 /* ------------------------------------------------------------------ */
4991 Unit *acc; /* -> accumulator Unit array */ in decMultiplyOp()
4993 void *allocacc=NULL; /* -> allocated accumulator, iff allocated */ in decMultiplyOp()
5006 #define FASTLAZY 18 /* carry resolution point [1->18] */ in decMultiplyOp()
5010 #define FASTLAZY 1844 /* carry resolution point [1->1844] */ in decMultiplyOp()
5016 uInt *zlhi=zlhibuff; /* -> lhs array */ in decMultiplyOp()
5017 uInt *alloclhi=NULL; /* -> allocated buffer, iff allocated */ in decMultiplyOp()
5019 uInt *zrhi=zrhibuff; /* -> rhs array */ in decMultiplyOp()
5020 uInt *allocrhi=NULL; /* -> allocated buffer, iff allocated */ in decMultiplyOp()
5023 uLong *zacc=zaccbuff; /* -> accumulator array for exact result */ in decMultiplyOp()
5027 uInt *lip, *rip; /* item pointers */ in decMultiplyOp()
5029 Int ilhs, irhs, iacc; /* item counts in the arrays */ in decMultiplyOp()
5041 decNumber *alloclhs=NULL; /* -> allocated buffer, iff allocated */ in decMultiplyOp()
5042 decNumber *allocrhs=NULL; /* -> allocated buffer, iff allocated */ in decMultiplyOp()
5050 bits=(uByte)((lhs->bits^rhs->bits)&DECNEG); in decMultiplyOp()
5058 if (((lhs->bits & DECINF)==0 && ISZERO(lhs)) in decMultiplyOp()
5059 ||((rhs->bits & DECINF)==0 && ISZERO(rhs))) { in decMultiplyOp()
5063 res->bits=bits|DECINF; /* infinity */ in decMultiplyOp()
5070 if (lhs->digits<rhs->digits) { /* swap... */ in decMultiplyOp()
5078 if (!set->extended) { in decMultiplyOp()
5080 if (lhs->digits>set->digits) { in decMultiplyOp()
5085 if (rhs->digits>set->digits) { in decMultiplyOp()
5098 if (rhs->digits>NEEDTWO) { /* use fastpath... */ in decMultiplyOp()
5100 ilhs=(lhs->digits+FASTDIGS-1)/FASTDIGS; /* [ceiling] */ in decMultiplyOp()
5101 irhs=(rhs->digits+FASTDIGS-1)/FASTDIGS; /* .. */ in decMultiplyOp()
5116 /* after the multiplication each 8-byte item becomes 9 1-byte */ in decMultiplyOp()
5133 acc=(Unit *)zacc; /* -> target Unit array */ in decMultiplyOp()
5139 for (count=lhs->digits, cup=lhs->lsu, lip=zlhi; count>0; lip++) in decMultiplyOp()
5141 p+=DECDPUN, cup++, count-=DECDPUN) in decMultiplyOp()
5143 lmsi=lip-1; /* save -> msi */ in decMultiplyOp()
5144 for (count=rhs->digits, cup=rhs->lsu, rip=zrhi; count>0; rip++) in decMultiplyOp()
5146 p+=DECDPUN, cup++, count-=DECDPUN) in decMultiplyOp()
5148 rmsi=rip-1; /* save -> msi */ in decMultiplyOp()
5156 /* Each uLong item in the accumulator can hold values up to */ in decMultiplyOp()
5157 /* 2**64-1, and each partial product can be as large as */ in decMultiplyOp()
5158 /* (10**FASTDIGS-1)**2. When FASTDIGS=9, this can be added to */ in decMultiplyOp()
5161 /* add -- every 162 digits. Similarly, when FASTDIGS=8, the */ in decMultiplyOp()
5171 for (rip=zrhi; rip<=rmsi; rip++) { /* over each item in rhs */ in decMultiplyOp()
5172 lp=zacc+(rip-zrhi); /* where to add the lhs */ in decMultiplyOp()
5173 for (lip=zlhi; lip<=lmsi; lip++, lp++) { /* over each item in lhs */ in decMultiplyOp()
5174 *lp+=(uLong)(*lip)*(*rip); /* [this should in-line] */ in decMultiplyOp()
5176 lazy--; in decMultiplyOp()
5183 /* lcarry can exceed 2**32-1, so check again; this check */ in decMultiplyOp()
5188 else { /* two-place carry [fairly rare] */ in decMultiplyOp()
5190 *(lp+2)+=carry2; /* add to item+2 */ in decMultiplyOp()
5191 *lp-=((uLong)FASTBASE*FASTBASE*carry2); /* [slow] */ in decMultiplyOp()
5192 carry=(uInt)(lcarry-((uLong)FASTBASE*carry2)); /* [inline] */ in decMultiplyOp()
5194 *(lp+1)+=carry; /* add to item above [inline] */ in decMultiplyOp()
5195 *lp-=((uLong)FASTBASE*carry); /* [inline] */ in decMultiplyOp()
5200 /* units. This can be done in-place in the accumulator and in */ in decMultiplyOp()
5201 /* 32-bit operations, because carries were resolved after the */ in decMultiplyOp()
5202 /* final add. This needs N-1 divides and multiplies for */ in decMultiplyOp()
5203 /* each item in the accumulator (which will become up to N */ in decMultiplyOp()
5206 uInt item=(uInt)*lp; /* decapitate to uInt */ in decMultiplyOp() local
5207 for (p=0; p<FASTDIGS-DECDPUN; p+=DECDPUN, up++) { in decMultiplyOp()
5208 uInt part=item/(DECDPUNMAX+1); in decMultiplyOp()
5209 *up=(Unit)(item-(part*(DECDPUNMAX+1))); in decMultiplyOp()
5210 item=part; in decMultiplyOp()
5212 *up=(Unit)item; up++; /* [final needs no division] */ in decMultiplyOp()
5214 accunits=up-acc; /* count of units */ in decMultiplyOp()
5220 acc=accbuff; /* -> assume buffer for accumulator */ in decMultiplyOp()
5221 needbytes=(D2U(lhs->digits)+D2U(rhs->digits))*sizeof(Unit); in decMultiplyOp()
5236 madlength=D2U(lhs->digits); /* this won't change */ in decMultiplyOp()
5237 mermsup=rhs->lsu+D2U(rhs->digits); /* -> msu+1 of multiplier */ in decMultiplyOp()
5239 for (mer=rhs->lsu; mer<mermsup; mer++) { in decMultiplyOp()
5241 /* If non-zero [optimization] add it... */ in decMultiplyOp()
5242 if (*mer!=0) accunits=decUnitAddSub(&acc[shift], accunits-shift, in decMultiplyOp()
5243 lhs->lsu, madlength, 0, in decMultiplyOp()
5256 /* common end-path */ in decMultiplyOp()
5264 res->bits=bits; /* set sign */ in decMultiplyOp()
5265 res->digits=decGetDigits(acc, accunits); /* count digits exactly */ in decMultiplyOp()
5267 /* There can be a 31-bit wrap in calculating the exponent. */ in decMultiplyOp()
5272 exponent=lhs->exponent+rhs->exponent; /* calculate exponent */ in decMultiplyOp()
5273 if (lhs->exponent<0 && rhs->exponent<0 && exponent>0) in decMultiplyOp()
5274 exponent=-2*DECNUMMAXE; /* force underflow */ in decMultiplyOp()
5275 res->exponent=exponent; /* OK to overwrite now */ in decMultiplyOp()
5279 decSetCoeff(res, set, acc, res->digits, &residue, status); in decMultiplyOp()
5295 /* ------------------------------------------------------------------ */
5296 /* decExpOp -- effect exponentiation */
5304 /* C must have space for set->digits digits. status is updated but */
5309 /* digits, emax, and -emin in the context must be less than */
5319 /* when A is a zero or -Infinity (giving 1 or 0 respectively). */
5320 /* ------------------------------------------------------------------ */
5325 /* pp79-91, ACM, June 1986. */
5329 /* extra variable-precision variables which are expensive in this */
5333 /* round-off error accumulation during the series evaluation. This */
5335 /* use Horner's scheme. Instead, the accumulation is done at double- */
5337 /* and do not accumulate round-off (and any round-off errors in the */
5350 /* exp(-x) where x can be the tiniest number (Ntiny). */
5358 /* sliding value (8-h), where potentially the range is reduced */
5362 /* limited by the cost of the raise-to-the power at the end, */
5365 /* x**10,000,000 needs 31 multiplications, all but one full-width. */
5369 /* 32-bit limits. */
5373 /* ------------------------------------------------------------------ */
5389 decNumber *allocrhs=NULL; /* non-NULL if rhs buffer allocated */ in decExpOp()
5391 /* the working precision will be no more than set->digits+8+1 */ in decExpOp()
5392 /* so for on-stack buffers DECBUFFER+9 is used, +1 in case DECBUFFER */ in decExpOp()
5397 decNumber *allocbuft=NULL; /* -> allocated buft, iff allocated */ in decExpOp()
5401 decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */ in decExpOp()
5417 if (decNumberIsNegative(rhs)) /* -Infinity -> +0 */ in decExpOp()
5419 else decNumberCopy(res, rhs); /* +Infinity -> self */ in decExpOp()
5424 if (ISZERO(rhs)) { /* zeros -> exact 1 */ in decExpOp()
5426 *res->lsu=1; /* .. */ in decExpOp()
5429 /* e**x when 0 < x < 0.66 is < 1+3x/2, hence can fast-path */ in decExpOp()
5431 /* 1. This also allows the later add-accumulate to always be */ in decExpOp()
5437 /* set->digits-1 zeros between the decimal point and the digit, */ in decExpOp()
5444 *d->lsu=4; /* set 4 .. */ in decExpOp()
5445 d->exponent=-set->digits; /* * 10**(-d) */ in decExpOp()
5446 if (decNumberIsNegative(rhs)) d->exponent--; /* negative case */ in decExpOp()
5452 Int shift=set->digits-1; in decExpOp()
5454 *res->lsu=1; /* .. */ in decExpOp()
5455 res->digits=decShiftToMost(res->lsu, 1, shift); in decExpOp()
5456 res->exponent=-shift; /* make 1.0000... */ in decExpOp()
5464 aset.emax=set->emax; /* usual bounds */ in decExpOp()
5465 aset.emin=set->emin; /* .. */ in decExpOp()
5470 h=rhs->exponent+rhs->digits; in decExpOp()
5474 /* overflow (or underflow to 0) is guaranteed -- so this case can */ in decExpOp()
5481 *a->lsu=2; /* not 1 but < exp(1) */ in decExpOp()
5482 if (decNumberIsNegative(rhs)) a->exponent=-2; /* make 0.02 */ in decExpOp()
5487 Int maxlever=(rhs->digits>8?1:0); in decExpOp()
5498 Int lever=MINI(8-h, maxlever); /* leverage attainable */ in decExpOp()
5499 Int use=-rhs->digits-lever; /* exponent to use for RHS */ in decExpOp()
5506 if (rhs->exponent!=use) { in decExpOp()
5508 needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit); in decExpOp()
5511 if (allocrhs==NULL) { /* hopeless -- abandon */ in decExpOp()
5517 newrhs->exponent=use; /* normalize; now <1 */ in decExpOp()
5528 /* this is set->digits+h+2. However, if x is 'over-precise' we */ in decExpOp()
5531 /* this case use x->digits+h+2 */ in decExpOp()
5532 p=MAXI(x->digits, set->digits)+h+2; /* [h<=8] */ in decExpOp()
5540 needbytes=sizeof(decNumber)+(D2U(p*2)-1)*sizeof(Unit); in decExpOp()
5543 if (allocbufa==NULL) { /* hopeless -- abandon */ in decExpOp()
5549 /* guaranteed to be larger than x->digits, so the initial copy */ in decExpOp()
5550 /* is safe); it may also be used for the raise-to-power */ in decExpOp()
5552 needbytes=sizeof(decNumber)+(D2U(p+2)-1)*sizeof(Unit); in decExpOp()
5555 if (allocbuft==NULL) { /* hopeless -- abandon */ in decExpOp()
5562 decNumberZero(a); *a->lsu=1; /* accumulator=1 */ in decExpOp()
5563 decNumberZero(d); *d->lsu=2; /* divisor=2 */ in decExpOp()
5590 if (((a->digits+a->exponent)>=(t->digits+t->exponent+p+1)) in decExpOp()
5591 && (a->digits>=p)) break; in decExpOp()
5599 iterations, *status, p, x->digits); in decExpOp()
5603 /* apply postconditioning: a=a**(10**h) -- this is calculated */ in decExpOp()
5606 Int seenbit=0; /* set once a 1-bit is seen */ in decExpOp()
5613 decNumberZero(t); *t->lsu=1; /* acc=1 */ in decExpOp()
5634 aset.digits=set->digits; /* [use default rounding] */ in decExpOp()
5646 /* ------------------------------------------------------------------ */
5647 /* Initial-estimate natural logarithm table */
5649 /* LNnn -- 90-entry 16-bit table for values from .10 through .99. */
5650 /* The result is a 4-digit encode of the coefficient (c=the */
5651 /* top 14 bits encoding 0-9999) and a 2-digit encode of the */
5652 /* exponent (e=the bottom 2 bits encoding 0-3) */
5656 /* v = -c * 10**(-e-3) */
5658 /* where e and c are extracted from entry k = LNnn[x-10] */
5661 /* ------------------------------------------------------------------ */
5674 /* ------------------------------------------------------------------ */
5675 /* decLnOp -- effect natural logarithm */
5683 /* C must have space for set->digits digits. */
5686 /* A<0 -> Invalid */
5687 /* A=0 -> -Infinity (Exact) */
5688 /* A=+Infinity -> +Infinity (Exact) */
5689 /* A=1 exactly -> 0 (Exact) */
5693 /* digits, emax, and -emin in the context must be less than */
5701 /* ------------------------------------------------------------------ */
5703 /* iteration calculating a' = a + x * exp(-a) - 1. See, for example, */
5706 /* The iteration ends when the adjustment x*exp(-a)-1 is tiny enough. */
5718 /* it typically needs 4-6 iterations for short numbers, and the */
5728 /* iterations will be needed -- so this would only speed up by */
5729 /* 20-25% and that probably does not justify increasing the table */
5734 /* ------------------------------------------------------------------ */
5748 decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */
5751 decNumber *allocbufb=NULL; /* -> allocated bufa, iff allocated */
5766 if (decNumberIsNegative(rhs)) /* -Infinity -> error */
5768 else decNumberCopy(res, rhs); /* +Infinity -> self */
5773 if (ISZERO(rhs)) { /* +/- zeros -> -Infinity */
5775 res->bits=DECINF|DECNEG; /* set - infinity */
5778 /* Non-zero negatives are bad... */
5779 if (decNumberIsNegative(rhs)) { /* -x -> error */
5786 if (rhs->exponent==0 && set->digits<=40) {
5788 if (rhs->lsu[0]==0 && rhs->lsu[1]==1 && rhs->digits==2) { /* ln(10) */
5790 if (rhs->lsu[0]==10 && rhs->digits==2) { /* ln(10) */
5797 if (rhs->lsu[0]==2 && rhs->digits==1) { /* ln(2) */
5807 /* the rhs is 'over-precise' then allow for all its digits to */
5809 /* digits are 9s) so in this case use rhs->digits+2. */
5810 p=MAXI(rhs->digits, MAXI(set->digits, 7))+2;
5812 /* Allocate space for the accumulator and the high-precision */
5815 /* rhs->digits+p digits. They are also made big enough for 16 */
5818 needbytes=sizeof(decNumber)+(D2U(MAXI(p,16))-1)*sizeof(Unit);
5821 if (allocbufa==NULL) { /* hopeless -- abandon */
5826 pp=p+rhs->digits;
5827 needbytes=sizeof(decNumber)+(D2U(MAXI(pp,16))-1)*sizeof(Unit);
5830 if (allocbufb==NULL) { /* hopeless -- abandon */
5845 decContextDefault(&aset, DEC_INIT_DECIMAL64); /* 16-digit extended */
5846 r=rhs->exponent+rhs->digits; /* 'normalised' exponent */
5849 b->exponent=-6; /* .. */
5856 b->exponent=0; /* make integer */
5858 if (t<10) t=X10(t); /* adjust single-digit b */
5859 t=LNnn[t-10]; /* look up ln(b) */
5861 b->exponent=-(t&3)-3; /* set exponent */
5862 b->bits=DECNEG; /* ln(0.10)->ln(0.99) always -ve */
5873 aset.emax=set->emax;
5874 aset.emin=set->emin;
5879 bset.emin=-DEC_MAX_MATH*2; /* adjustment calculation */
5885 /* 34+2, which is ideal for standard-sized numbers] */
5887 bset.digits=pp+rhs->digits; /* wider context */
5893 /* calculate the adjustment (exp(-a)*x-1) into b. This is a */
5897 /* range for calculating exp(-a) when a is the tiniest subnormal. */
5898 a->bits^=DECNEG; /* make -a */
5899 decExpOp(b, a, &bset, &ignore); /* b=exp(-a) */
5900 a->bits^=DECNEG; /* restore sign of a */
5903 decAddOp(b, b, &numone, &bset, DECNEG, &ignore); /* b=b-1 */
5908 /* set->digits+1 digits (or it is zero) -- this is a looser */
5914 (a->digits+a->exponent)>=(b->digits+b->exponent+set->digits+1)) {
5915 if (a->digits==p) break;
5918 if (cmp.lsu[0]==0) a->exponent=0; /* yes, exact 0 */
5923 if (decNumberIsZero(b)) b->exponent=a->exponent-p;
5933 bset.digits=pp+rhs->digits; /* wider context */
5937 /* just a sanity check; remove the test to show always */
5940 iterations, *status, p, rhs->digits);
5946 aset.digits=set->digits; /* [use default rounding] */
5957 /* ------------------------------------------------------------------ */
5958 /* decQuantizeOp -- force exponent to requested value */
5961 /* of C (by rounding or shifting) such that the exponent (-scale) */
5974 /* C must have space for set->digits digits. */
5978 /* ------------------------------------------------------------------ */
5983 decNumber *alloclhs=NULL; /* non-NULL if rounded lhs allocated */
5987 Int reqdigits=set->digits; /* requested DIGITS */
5988 Int reqexp; /* requested exponent [-scale] */
5990 Int etiny=set->emin-(reqdigits-1);
5998 if (!set->extended) {
6000 if (lhs->digits>reqdigits) {
6005 if (rhs->digits>reqdigits) { /* [this only checks lostDigits] */
6020 else if ((lhs->bits ^ rhs->bits) & DECINF)
6028 if (quant) reqexp=inrhs->exponent; /* quantize -- match exponents */
6029 else { /* rescale -- use value of rhs */
6031 /* which could be from -1999999997 to +999999999, thanks to */
6037 if (!set->extended) etiny=set->emin; /* no subnormals */
6043 || (reqexp>set->emax)) { /* > emax */
6050 res->exponent=reqexp; /* .. just set exponent */
6052 if (!set->extended) res->bits=0; /* subset specification; no -0 */
6055 else { /* non-zero lhs */
6056 Int adjust=reqexp-lhs->exponent; /* digit adjustment needed */
6058 if ((lhs->digits-adjust)>reqdigits) {
6068 workset.digits=lhs->digits-adjust; /* set requested length */
6075 if (res->exponent>reqexp) {
6076 /* re-check needed, e.g., for quantize(0.9999, 0.001) under */
6077 /* set->digits==3 */
6078 if (res->digits==reqdigits) { /* cannot shift by 1 */
6083 res->digits=decShiftToMost(res->lsu, res->digits, 1); /* shift */
6084 res->exponent--; /* (re)adjust the exponent. */
6087 if (ISZERO(res) && !set->extended) res->bits=0; /* subset; no -0 */
6091 /* this will increase the length of the coefficient by -adjust */
6097 res->digits=decShiftToMost(res->lsu, res->digits, -adjust);
6098 res->exponent+=adjust; /* adjust the exponent */
6101 } /* non-zero */
6105 if (res->exponent>set->emax-res->digits+1) { /* too big */
6122 /* ------------------------------------------------------------------ */
6123 /* decCompareOp -- compare, min, or max two Numbers */
6126 /* COMPARE -- returns the signum (as a number) giving the */
6129 /* COMPSIG -- as COMPARE except that a quiet NaN raises */
6131 /* COMPMAX -- returns the larger of the operands, using the */
6133 /* COMPMAXMAG -- ditto, comparing absolute values */
6134 /* COMPMIN -- the 754r minnum operation */
6135 /* COMPMINMAG -- ditto, comparing absolute values */
6136 /* COMTOTAL -- returns the signum (as a number) giving the */
6146 /* C must have space for one digit for COMPARE or set->digits for */
6148 /* ------------------------------------------------------------------ */
6151 /* ------------------------------------------------------------------ */
6156 decNumber *alloclhs=NULL; /* non-NULL if rounded lhs allocated */
6168 if (!set->extended) {
6170 if (lhs->digits>set->digits) {
6175 if (rhs->digits>set->digits) {
6187 result=-1;
6198 merged=(lhs->bits | rhs->bits) & (DECSNAN | DECNAN);
6206 if (!decNumberIsNaN(lhs)) result=-1;
6209 else if (decNumberIsSNaN(lhs) && decNumberIsQNaN(rhs)) result=-1;
6213 result=decUnitCompare(lhs->lsu, D2U(lhs->digits),
6214 rhs->lsu, D2U(rhs->digits), 0);
6217 if (decNumberIsNegative(lhs)) result=-result;
6221 else if (merged & DECSNAN); /* sNaN -> qNaN */
6223 /* min or max -- 754r rules ignore single NaN */
6225 /* just one NaN; force choice to be the non-NaN operand */
6227 if (lhs->bits & DECNAN) result=-1; /* pick rhs */
6247 if (lhs->exponent!=rhs->exponent) {
6248 if (lhs->exponent<rhs->exponent) result=-1;
6250 if (decNumberIsNegative(lhs)) result=-result;
6252 } /* total-order by exponent */
6254 if (result!=0) { /* must be -1 or +1 */
6255 *res->lsu=1;
6256 if (result<0) res->bits=DECNEG;
6260 else { /* MAX or MIN, non-NaN result */
6266 uByte slhs=(lhs->bits & DECNEG);
6267 uByte srhs=(rhs->bits & DECNEG);
6269 if (!set->extended) { /* subset: force left-hand */
6276 if (slhs) result=-1; /* rhs is max */
6280 if (lhs->exponent<rhs->exponent) result=+1;
6281 else result=-1;
6285 if (lhs->exponent>rhs->exponent) result=+1;
6286 else result=-1;
6290 /* here result will be non-0; reverse if looking for MIN */
6291 if (op==COMPMIN || op==COMPMINMAG) result=-result;
6305 /* ------------------------------------------------------------------ */
6306 /* decCompare -- compare two decNumbers by numerical value */
6312 /* Arg3 is 1 for a sign-independent compare, 0 otherwise */
6314 /* returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure */
6316 /* ------------------------------------------------------------------ */
6327 /* RHS is non-zero */
6328 if (result==0) return -1; /* LHS is 0; RHS wins */
6329 /* [here, both non-zero, result=1] */
6332 if (result && decNumberIsNegative(lhs)) result=-1;
6335 else if (decNumberIsNegative(rhs)) sigr=-1;
6337 if (result < sigr) return -1; /* L < R, return -1 */
6341 /* signums are the same; both are non-zero */
6342 if ((lhs->bits | rhs->bits) & DECINF) { /* one or more infinities */
6345 else result=-result; /* only rhs infinite */
6350 if (lhs->exponent>rhs->exponent) { /* LHS exponent larger */
6355 result=-result;
6357 compare=decUnitCompare(lhs->lsu, D2U(lhs->digits),
6358 rhs->lsu, D2U(rhs->digits),
6359 rhs->exponent-lhs->exponent);
6364 /* ------------------------------------------------------------------ */
6365 /* decUnitCompare -- compare two >=0 integers in Unit arrays */
6369 /* B has an exponent of E (which must be non-negative) */
6377 /* returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure */
6380 /* ------------------------------------------------------------------ */
6385 Unit *allocacc=NULL; /* -> allocated acc buffer, iff allocated */
6392 if (alength<blength) return -1;
6393 /* same number of units in both -- need unit-by-unit compare */
6394 l=a+alength-1;
6395 r=b+alength-1;
6396 for (;l>=a; l--, r--) {
6398 if (*l<*r) return -1;
6406 if (alength+1<blength+(Int)D2U(exp)) return -1;
6417 if (allocacc==NULL) return BADINT; /* hopeless -- abandon */
6423 /* subtract [A+B*(-m)] */
6425 -(Int)powers[exprem]);
6427 if (accunits<0) result=-1; /* negative result */
6428 else { /* non-negative result */
6430 for (u=acc; u<acc+accunits-1 && *u==0;) u++;
6438 /* ------------------------------------------------------------------ */
6439 /* decUnitAddSub -- add or subtract two >=0 integers in Unit arrays */
6445 /* Where M is in the range -DECDPUNMAX through +DECDPUNMAX. */
6467 /* returns the count of Units written to C, which will be non-zero */
6473 /* safe, allowing space if necessary for a one-Unit carry. */
6475 /* This routine is severely performance-critical; *any* change here */
6477 /* In particular, trickery here tends to be counter-productive, as */
6479 /* register-poor architectures. Avoiding divisions is nearly */
6484 /* ------------------------------------------------------------------ */
6530 carry+=((eInt)*b)*m; /* [special-casing m=1/-1 */
6532 /* here carry is new Unit of digits; it could be +ve or -ve */
6533 if ((ueInt)carry<=DECDPUNMAX) { /* fastpath 0-DECDPUNMAX */
6538 #if DECDPUN==4 /* use divide-by-multiply */
6541 *c=(Unit)(carry-est*(DECDPUNMAX+1)); /* remainder */
6545 *c-=DECDPUNMAX+1;
6551 *c=(Unit)(carry-est*(DECDPUNMAX+1));
6552 carry=est-(DECDPUNMAX+1); /* correctly negative */
6555 *c-=DECDPUNMAX+1;
6559 *c=(Unit)(carry-est*(DECDPUNMAX+1)); /* remainder */
6563 *c-=DECDPUNMAX+1;
6569 *c=(Unit)(carry-est*(DECDPUNMAX+1));
6570 carry=est-(DECDPUNMAX+1); /* correctly negative */
6573 *c-=DECDPUNMAX+1;
6578 *c=(Unit)(carry-est*(DECDPUNMAX+1)); /* remainder */
6585 *c=(Unit)(carry-est*(DECDPUNMAX+1));
6586 carry=est-(DECDPUNMAX+1); /* correctly negative */
6590 *c=(Unit)(carry-(DECDPUNMAX+1)); /* [helps additions] */
6602 carry=carry/(DECDPUNMAX+1)-(DECDPUNMAX+1);
6617 /* here carry is new Unit of digits; it could be +ve or -ve and */
6619 if ((ueInt)carry<=DECDPUNMAX) { /* fastpath 0-DECDPUNMAX */
6625 #if DECDPUN==4 /* use divide-by-multiply */
6628 *c=(Unit)(carry-est*(DECDPUNMAX+1)); /* remainder */
6632 *c-=DECDPUNMAX+1;
6638 *c=(Unit)(carry-est*(DECDPUNMAX+1));
6639 carry=est-(DECDPUNMAX+1); /* correctly negative */
6642 *c-=DECDPUNMAX+1;
6646 *c=(Unit)(carry-est*(DECDPUNMAX+1)); /* remainder */
6650 *c-=DECDPUNMAX+1;
6656 *c=(Unit)(carry-est*(DECDPUNMAX+1));
6657 carry=est-(DECDPUNMAX+1); /* correctly negative */
6660 *c-=DECDPUNMAX+1;
6664 *c=(Unit)(carry-est*(DECDPUNMAX+1)); /* remainder */
6671 *c=(Unit)(carry-est*(DECDPUNMAX+1));
6672 carry=est-(DECDPUNMAX+1); /* correctly negative */
6675 *c=(Unit)(carry-(DECDPUNMAX+1));
6688 carry=carry/(DECDPUNMAX+1)-(DECDPUNMAX+1);
6694 if (carry==0) return c-clsu; /* no carry, so no more to do */
6698 return c-clsu;
6700 /* -ve carry: it's a borrow; complement needed */
6703 add=DECDPUNMAX+add-*c;
6713 /* add an extra unit iff it would be non-zero */
6717 if ((add-carry-1)!=0) {
6718 *c=(Unit)(add-carry-1);
6721 return clsu-c; /* -ve result indicates borrowed */
6724 /* ------------------------------------------------------------------ */
6725 /* decTrim -- trim trailing zeros or normalize */
6729 /* all is 1 to remove all trailing zeros, 0 for just fraction ones */
6737 /* ------------------------------------------------------------------ */
6742 Unit *up; /* -> current Unit */
6749 if ((dn->bits & DECSPECIAL) /* fast exit if special .. */
6750 || (*dn->lsu & 0x01)) return dn; /* .. or odd */
6752 dn->exponent=0; /* (sign is preserved) */
6757 exp=dn->exponent;
6758 cut=1; /* digit (1-DECDPUN) in Unit */
6759 up=dn->lsu; /* -> current Unit */
6760 for (d=0; d<dn->digits-1; d++) { /* [don't strip the final digit] */
6764 if ((*up-quot*powers[cut])!=0) break; /* found non-0 digit */
6766 if (*up%powers[cut]!=0) break; /* found non-0 digit */
6785 if (set->clamp) {
6786 Int maxd=set->emax-set->digits+1-dn->exponent;
6792 decShiftToLeast(dn->lsu, D2U(dn->digits), d);
6793 dn->exponent+=d; /* maintain numerical value */
6794 dn->digits-=d; /* new length */
6799 /* ------------------------------------------------------------------ */
6800 /* decReverse -- reverse a Unit array in place */
6808 /* ------------------------------------------------------------------ */
6811 for (; ulo<uhi; ulo++, uhi--) {
6819 /* ------------------------------------------------------------------ */
6820 /* decShiftToMost -- shift digits in array towards most significant */
6831 /* ------------------------------------------------------------------ */
6838 if ((digits+shift)<=DECDPUN) { /* [fastpath] single-unit case */
6844 source=uar+D2U(digits)-1; /* where msu comes from */
6846 cut=DECDPUN-MSUDIGITS(shift); /* where to slice */
6847 if (cut==0) { /* unit-boundary case */
6848 for (; source>=uar; source--, target--) *target=*source;
6851 first=uar+D2U(digits+shift)-1; /* where msu of source will end up */
6852 for (; source>=uar; source--, target--) {
6856 uInt rem=*source-quot*powers[cut];
6863 next=rem*powers[DECDPUN-cut]; /* save remainder for next Unit */
6865 } /* shift-move */
6868 for (; target>=uar; target--) {
6875 /* ------------------------------------------------------------------ */
6876 /* decShiftToLeast -- shift digits in array towards least significant */
6880 /* shift is the number of digits to remove from the lsu end; it */
6887 /* ------------------------------------------------------------------ */
6901 if (cut==DECDPUN) { /* unit-boundary case; easy */
6904 return target-uar;
6908 up=uar+D2U(shift-cut); /* source; correct to whole Units */
6909 count=units*DECDPUN-shift; /* the maximum new length */
6917 count-=(DECDPUN-cut);
6923 rem=*up-quot*powers[cut];
6928 *target=(Unit)(*target+rem*powers[DECDPUN-cut]);
6929 count-=cut;
6932 return target-uar+1;
6936 /* ------------------------------------------------------------------ */
6937 /* decRoundOperand -- round an operand [used for subset only] */
6939 /* dn is the number to round (dn->digits is > set->digits) */
6953 /* ------------------------------------------------------------------ */
6963 +(D2U(set->digits)-1)*sizeof(Unit));
6978 /* ------------------------------------------------------------------ */
6979 /* decCopyFit -- copy a number, truncating the coefficient if needed */
6987 /* (dest==src is allowed and will be a no-op if fits) */
6989 /* ------------------------------------------------------------------ */
6992 dest->bits=src->bits;
6993 dest->exponent=src->exponent;
6994 decSetCoeff(dest, set, src->lsu, src->digits, residue, status);
6997 /* ------------------------------------------------------------------ */
6998 /* decSetCoeff -- set the coefficient of a number */
7001 /* It must have space for set->digits digits */
7003 /* lsu -> lsu of the source coefficient [may be dn->lsu] */
7004 /* len is digits in the source coefficient [may be dn->digits] */
7014 /* dn->lsu and len must == dn->digits. */
7016 /* Note that the coefficient length (len) may be < set->digits, and */
7017 /* in this case this merely copies the coefficient (or is a no-op */
7018 /* if dn->lsu==lsu). */
7021 /* decSetSubnormal) the value of set->digits may be less than one, */
7025 /* dn->digits, dn->lsu (and as required), and dn->exponent are */
7026 /* updated as necessary. dn->bits (sign) is unchanged. */
7029 /* DEC_Inexact status is set if any non-zero digits are discarded, or */
7030 /* incoming residue was non-0 (implies rounded) */
7031 /* ------------------------------------------------------------------ */
7032 /* mapping array: maps 0-9 to canonical residues, so that a residue */
7033 /* can be adjusted in the range [-1, +1] and achieve correct rounding */
7047 discard=len-set->digits; /* digits to discard */
7049 if (dn->lsu!=lsu) { /* copy needed */
7053 for (target=dn->lsu; count>0; target++, up++, count-=DECDPUN)
7055 dn->digits=len; /* set the new length */
7057 /* dn->exponent and residue are unchanged, record any inexactitude */
7063 dn->exponent+=discard; /* maintain numerical value */
7072 for (up=lsu; count>0; up++, count-=DECDPUN) if (*up!=0) { /* found non-0 */
7078 *dn->lsu=0; /* coefficient will now be 0 */
7079 dn->digits=1; /* .. */
7096 /* here up -> Unit with first discarded digit */
7097 cut=discard-(count-DECDPUN)-1;
7098 if (cut==DECDPUN-1) { /* unit-boundary case (fast) */
7108 if (set->digits<=0) { /* special for Quantize/Subnormal :-( */
7109 *dn->lsu=0; /* .. result is 0 */
7110 dn->digits=1; /* .. */
7113 count=set->digits; /* now digits to end up with */
7114 dn->digits=count; /* set the new length */
7116 /* on unit boundary, so shift-down copy loop is simple */
7117 for (target=dn->lsu; count>0; target++, up++, count-=DECDPUN)
7120 } /* unit-boundary case */
7129 rem=*up-quot*powers[cut];
7140 discard1=quot-X10(temp);
7150 /* here: up -> Unit of the array with bottom digit */
7152 /* quot holds the uncut high-order digits for the current unit */
7153 if (set->digits<=0) { /* special for Quantize/Subnormal :-( */
7154 *dn->lsu=0; /* .. result is 0 */
7155 dn->digits=1; /* .. */
7158 count=set->digits; /* now digits to end up with */
7159 dn->digits=count; /* set the new length */
7160 /* shift-copy the coefficient array to the result number */
7161 for (target=dn->lsu; ; target++) {
7163 count-=(DECDPUN-cut);
7169 rem=*up-quot*powers[cut];
7174 *target=(Unit)(*target+rem*powers[DECDPUN-cut]);
7175 count-=cut;
7177 } /* shift-copy loop */
7185 /* ------------------------------------------------------------------ */
7186 /* decApplyRound -- apply pending rounding to a number */
7188 /* dn is the number, with space for set->digits digits */
7192 /* 6-9: rounding digit is >5 */
7193 /* 5: rounding digit is exactly half-way */
7194 /* 1-4: rounding digit is <5 and >0 */
7196 /* -1: as 1, but the hidden digits are subtractive, that */
7198 /* coefficient must be non-0. This case occurs when */
7207 /* -- the coefficient was increased and is all nines (in which */
7209 /* the caller does not need to re-test for overflow) */
7211 /* -- the coefficient was decreased and becomes all nines (in which */
7216 /* ------------------------------------------------------------------ */
7220 /* -1 if coefficient needs to be decremented */
7227 switch (set->round) {
7233 /* case it is bumped down and then up -- a no-op) */
7234 Int lsd5=*dn->lsu%5; /* get lsd and quintate */
7235 if (residue<0 && lsd5!=1) bump=-1;
7238 break;} /* r-05 */
7242 if (residue<0) bump=-1;
7243 break;} /* r-d */
7247 break;} /* r-h-d */
7253 if (*dn->lsu & 0x01) bump=1;
7255 break;} /* r-h-e */
7259 break;} /* r-h-u */
7263 break;} /* r-u */
7269 if (residue<0) bump=-1;
7274 break;} /* r-c */
7280 if (residue<0) bump=-1;
7285 break;} /* r-f */
7290 printf("Unknown rounding mode: %d\n", set->round);
7302 /* Similarly handle all-nines result if bumping down. */
7305 uInt count=dn->digits; /* digits to be checked */
7306 for (up=dn->lsu; ; up++) {
7309 if (*up!=powers[count]-1) break; /* not still 9s */
7311 *up=(Unit)powers[count-1]; /* here 999 -> 100 etc. */
7312 for (up=up-1; up>=dn->lsu; up--) *up=0; /* others all to 0 */
7313 dn->exponent++; /* and bump exponent */
7315 if ((dn->exponent+dn->digits)>set->emax+1) {
7322 count-=DECDPUN;
7325 else { /* -1 */
7326 /* here checking for a pre-bump of 1000... (leading 1, all */
7329 uInt count=dn->digits; /* digits to be checked */
7330 for (up=dn->lsu; ; up++) {
7333 if (*up!=powers[count-1]) break; /* not 100.. */
7336 *up=(Unit)powers[count]-1; /* here 100 in msu -> 999 */
7337 /* others all to all-nines, too */
7338 for (up=up-1; up>=dn->lsu; up--) *up=(Unit)powers[DECDPUN]-1;
7339 dn->exponent--; /* and bump exponent */
7344 /* printf(">> emin=%d exp=%d sdig=%d\n", set->emin, */
7345 /* dn->exponent, set->digits); */
7346 if (dn->exponent+1==set->emin-set->digits+1) {
7347 if (count==1 && dn->digits==1) *sup=0; /* here 9 -> 0[.9] */
7349 *sup=(Unit)powers[count-1]-1; /* here 999.. in msu -> 99.. */
7350 dn->digits--;
7352 dn->exponent++;
7360 count-=DECDPUN;
7366 decUnitAddSub(dn->lsu, D2U(dn->digits), uarrone, 1, 0, dn->lsu, bump);
7370 /* ------------------------------------------------------------------ */
7371 /* decFinish -- finish processing a number */
7385 /* ------------------------------------------------------------------ */
7388 if (!set->extended) {
7390 dn->exponent=0; /* clean exponent .. */
7391 dn->bits=0; /* .. and sign */
7394 if (dn->exponent>=0) { /* non-negative exponent */
7396 if (set->digits >= (dn->exponent+dn->digits)) {
7397 dn->digits=decShiftToMost(dn->lsu, dn->digits, dn->exponent);
7398 dn->exponent=0;
7407 /* ------------------------------------------------------------------ */
7408 /* decFinalize -- final check, clamp, and round of a number */
7420 /* ------------------------------------------------------------------ */
7424 Int tinyexp=set->emin-dn->digits+1; /* precalculate subnormal boundary */
7432 if (dn->exponent<=tinyexp) { /* prefilter */
7436 if (dn->exponent<tinyexp) {
7444 nmin.exponent=set->emin;
7461 if (dn->exponent<=set->emax-set->digits+1) return; /* neither needed */
7465 if (dn->exponent>set->emax-dn->digits+1) { /* too big */
7470 if (!set->clamp) return;
7472 /* here when need to apply the IEEE exponent clamp (fold-down) */
7473 shift=dn->exponent-(set->emax-set->digits+1);
7475 /* shift coefficient (if non-zero) */
7477 dn->digits=decShiftToMost(dn->lsu, dn->digits, shift);
7479 dn->exponent-=shift; /* adjust the exponent to match */
7484 /* ------------------------------------------------------------------ */
7485 /* decSetOverflow -- set number to proper overflow value */
7494 /* ------------------------------------------------------------------ */
7497 uByte sign=dn->bits&DECNEG; /* clean and save sign bit */
7500 Int emax=set->emax; /* limit value */
7501 if (set->clamp) emax-=set->digits-1; /* lower if clamping */
7502 if (dn->exponent>emax) { /* clamp required */
7503 dn->exponent=emax;
7510 switch (set->round) {
7513 break;} /* r-d */
7516 break;} /* r-05 */
7518 if (sign) needmax=1; /* Infinity if non-negative */
7519 break;} /* r-c */
7522 break;} /* r-f */
7527 dn->bits=sign; /* set sign */
7529 else dn->bits=sign|DECINF; /* Value is +/-Infinity */
7533 /* ------------------------------------------------------------------ */
7534 /* decSetMaxValue -- set number to +Nmax (maximum normal value) */
7540 /* ------------------------------------------------------------------ */
7543 Int count=set->digits; /* nines to add */
7544 dn->digits=count;
7546 for (up=dn->lsu; ; up++) {
7549 *up=(Unit)(powers[count]-1);
7552 count-=DECDPUN; /* filled those digits */
7554 dn->bits=0; /* + sign */
7555 dn->exponent=set->emax-set->digits+1;
7558 /* ------------------------------------------------------------------ */
7559 /* decSetSubnormal -- process value whose exponent is <Emin */
7574 /* ------------------------------------------------------------------ */
7582 if (!set->extended) {
7590 /* Full arithmetic -- allow subnormals, rounded to minimum exponent */
7592 etiny=set->emin-(set->digits-1); /* smallest allowed exponent */
7595 /* residue can never be non-zero here */
7602 if (dn->exponent<etiny) { /* clamp required */
7603 dn->exponent=etiny;
7609 *status|=DEC_Subnormal; /* have a non-zero subnormal */
7610 adjust=etiny-dn->exponent; /* calculate digits to remove */
7612 /* residue can never be non-zero here, except in the Nmin-residue */
7613 /* case (which is a subnormal result), so can take fast-path here */
7622 workset.digits=dn->digits-adjust; /* set requested length */
7623 workset.emin-=adjust; /* and adjust emin to match */
7625 decSetCoeff(dn, &workset, dn->lsu, dn->digits, residue, status);
7634 if (dn->exponent>etiny) {
7635 dn->digits=decShiftToMost(dn->lsu, dn->digits, 1);
7636 dn->exponent--; /* (re)adjust the exponent. */
7643 /* ------------------------------------------------------------------ */
7644 /* decCheckMath - check entry conditions for a math function */
7652 /* returns non-zero if status is changed, 0 otherwise */
7656 /* digits, emax, and -emin in the context must be less than */
7658 /* non-zero. Invalid_operation is set in the status if a */
7660 /* ------------------------------------------------------------------ */
7664 if (set->digits>DEC_MAX_MATH
7665 || set->emax>DEC_MAX_MATH
7666 || -set->emin>DEC_MAX_MATH) *status|=DEC_Invalid_context;
7667 else if ((rhs->digits>DEC_MAX_MATH
7668 || rhs->exponent+rhs->digits>DEC_MAX_MATH+1
7669 || rhs->exponent+rhs->digits<2*(1-DEC_MAX_MATH))
7674 /* ------------------------------------------------------------------ */
7675 /* decGetInt -- get integer from a number */
7680 /* BADINT if there is a non-zero fraction */
7688 /* ------------------------------------------------------------------ */
7693 Int ilength=dn->digits+dn->exponent; /* integral length */
7694 Flag neg=decNumberIsNegative(dn); /* 1 if -ve */
7701 #if DEC_MIN_EMIN < -999999999
7706 up=dn->lsu; /* ready for lsu */
7708 if (dn->exponent>=0) { /* relatively easy */
7710 got=dn->exponent;
7712 else { /* -ve exponent; some fractional part to check and discard */
7713 Int count=-dn->exponent; /* digits to discard */
7716 if (*up!=0) return BADINT; /* non-zero Unit to discard */
7717 count-=DECDPUN;
7722 /* slice off fraction digits and check for non-zero */
7725 rem=*up-theInt*powers[count];
7730 if (rem!=0) return BADINT; /* non-zero fraction */
7732 got=DECDPUN-count; /* number of digits so far */
7749 if (theInt/(Int)powers[got-DECDPUN]!=(Int)*(up-1)) ilength=11;
7762 if (neg) theInt=-theInt; /* apply sign */
7766 /* ------------------------------------------------------------------ */
7767 /* decDecap -- decapitate the coefficient of a number */
7771 /* this must be <= dn->digits (if equal, the coefficient is */
7774 /* Returns dn; dn->digits will be <= the initial digits less drop */
7776 /* which will also be removed). Only dn->lsu and dn->digits change. */
7777 /* ------------------------------------------------------------------ */
7779 Unit *msu; /* -> target cut point */
7781 if (drop>=dn->digits) { /* losing the whole thing */
7783 if (drop>dn->digits)
7785 (LI)drop, (LI)dn->digits);
7787 dn->lsu[0]=0;
7788 dn->digits=1;
7791 msu=dn->lsu+D2U(dn->digits-drop)-1; /* -> likely msu */
7792 cut=MSUDIGITS(dn->digits-drop); /* digits to be in use in msu */
7795 dn->digits=decGetDigits(dn->lsu, msu-dn->lsu+1);
7799 /* ------------------------------------------------------------------ */
7800 /* decBiStr -- compare string with pairwise options */
7813 /* ------------------------------------------------------------------ */
7823 /* ------------------------------------------------------------------ */
7824 /* decNaNs -- handle NaN operand or operands */
7836 /* ------------------------------------------------------------------ */
7842 if (lhs->bits & DECSNAN)
7845 else if (rhs->bits & DECSNAN) {
7849 else if (lhs->bits & DECNAN);
7853 if (lhs->digits<=set->digits) decNumberCopy(res, lhs); /* easy */
7858 res->bits=lhs->bits; /* need sign etc. */
7859 uresp1=res->lsu+D2U(set->digits);
7860 for (ur=res->lsu, ul=lhs->lsu; ur<uresp1; ur++, ul++) *ur=*ul;
7861 res->digits=D2U(set->digits)*DECDPUN;
7863 if (res->digits>set->digits) decDecap(res, res->digits-set->digits);
7866 res->bits&=~DECSNAN; /* convert any sNaN to NaN, while */
7867 res->bits|=DECNAN; /* .. preserving sign */
7868 res->exponent=0; /* clean exponent */
7873 /* ------------------------------------------------------------------ */
7874 /* decStatus -- apply non-zero status */
7881 /* unless the error was an overflow, divide-by-zero, or underflow, */
7887 /* ------------------------------------------------------------------ */
7889 if (status & DEC_NaNs) { /* error status -> NaN */
7894 dn->bits=DECNAN; /* and make a quiet NaN */
7901 /* ------------------------------------------------------------------ */
7902 /* decGetDigits -- count digits in a Units array */
7912 /* ------------------------------------------------------------------ */
7915 Unit *up=uar+(len-1); /* -> msu */
7916 Int digits=(len-1)*DECDPUN+1; /* possible digits excluding msu */
7925 for (; up>=uar; up--) {
7928 digits-=DECDPUN; /* adjust for 0 unit */
7930 /* found the first (most significant) non-zero Unit */
7932 if (*up<10) break; /* is 1-9 */
7935 if (*up<100) break; /* is 10-99 */
7938 if (*up<1000) break; /* is 100-999 */
7951 /* ------------------------------------------------------------------ */
7952 /* mulUInt128ByPowOf10 -- multiply a 128-bit unsigned integer by a */
7955 /* The 128-bit factor composed of plow and phigh is multiplied */
7965 /* ------------------------------------------------------------------ */
7969 if (mulu128(plow, phigh, powers[ARRAY_SIZE(powers) - 1])) {
7973 pow10 -= ARRAY_SIZE(powers) - 1;
7984 /* ------------------------------------------------------------------ */
7985 /* decNumberShow -- display a number [debug aid] */
7989 /* or: sign, special-value */
7990 /* ------------------------------------------------------------------ */
8000 if (decNumberIsNegative(dn)) isign='-';
8002 if (dn->bits&DECSPECIAL) { /* Is a special value */
8005 if (dn->bits&DECSNAN) printf("sNaN"); /* signalling NaN */
8009 if (dn->exponent==0 && dn->digits==1 && *dn->lsu==0) {
8017 up=dn->lsu+D2U(dn->digits)-1; /* msu */
8019 for (up=up-1; up>=dn->lsu; up--) {
8022 for (cut=DECDPUN-1; cut>=0; cut--) {
8024 u-=d*powers[cut];
8028 if (dn->exponent!=0) {
8030 if (dn->exponent<0) esign='-';
8031 printf(" E%c%ld", esign, (LI)abs(dn->exponent));
8033 printf(" [%ld]\n", (LI)dn->digits);
8038 /* ------------------------------------------------------------------ */
8039 /* decDumpAr -- display a unit array [debug/check aid] */
8040 /* name is a single-character tag name */
8043 /* ------------------------------------------------------------------ */
8067 for (i=len-1; i>=0; i--) {
8068 if (i==len-1) printf("%ld ", (LI)ar[i]);
8076 /* ------------------------------------------------------------------ */
8077 /* decCheckOperands -- check operand(s) to a routine */
8088 /* ------------------------------------------------------------------ */
8099 && (set->digits<1 || set->round>=DEC_ROUND_MAX)) {
8103 (LI)set->digits, (LI)set->round);
8121 res->bits=DECNAN; /* qNaN */
8127 /* ------------------------------------------------------------------ */
8128 /* decCheckNumber -- check a number */
8134 /* ------------------------------------------------------------------ */
8149 if (dn->bits & DECSPECIAL) {
8150 if (dn->exponent!=0) {
8153 (LI)dn->exponent, dn->bits);
8159 if (dn->digits!=1) {
8161 printf("Digits %ld (not 1) for an infinity.\n", (LI)dn->digits);
8164 if (*dn->lsu!=0) {
8166 printf("LSU %ld (not 0) for an infinity.\n", (LI)*dn->lsu);
8168 decDumpAr('I', dn->lsu, D2U(dn->digits));
8177 if (dn->digits<1 || dn->digits>DECNUMMAXP) {
8179 printf("Digits %ld in number.\n", (LI)dn->digits);
8183 d=dn->digits;
8185 for (up=dn->lsu; d>0; up++) {
8188 maxuint=powers[d]-1;
8189 if (dn->digits>1 && *up<powers[d-1]) {
8198 printf("Bad Unit [%08lx] in %ld-digit number at offset %ld [maxuint %ld].\n",
8199 (LI)*up, (LI)dn->digits, (LI)(up-dn->lsu), (LI)maxuint);
8202 d-=DECDPUN;
8206 /* which are out of the set->emin/set->emax and set->digits range */
8207 /* (just as they can have more digits than set->digits). */
8208 ae=dn->exponent+dn->digits-1; /* adjusted exponent */
8212 if (ae<emin-(digits-1)) {
8228 /* ------------------------------------------------------------------ */
8229 /* decCheckInexact -- check a normal finite inexact result has digits */
8235 /* ------------------------------------------------------------------ */
8238 if ((set->status & (DEC_Inexact|DEC_Subnormal))==DEC_Inexact
8239 && (set->digits!=dn->digits) && !(dn->bits & DECSPECIAL)) {
8242 (LI)dn->digits);
8249 if (dn!=NULL && dn->digits==0) set->status|=DEC_Invalid_operation;
8258 /* ------------------------------------------------------------------ */
8259 /* decMalloc -- accountable allocation routine */
8265 /* ------------------------------------------------------------------ */
8268 /* 0-3 the original length requested */
8269 /* 4-7 buffer corruption detection fence (DECFENCE, x4) */
8271 /* ------------------------------------------------------------------ */
8274 void *alloc; /* -> allocated storage */
8278 alloc=malloc(size); /* -> allocated storage */
8282 j=(uInt *)alloc; /* -> first four bytes */
8287 return b0+8; /* -> play area */
8290 /* ------------------------------------------------------------------ */
8291 /* decFree -- accountable free routine */
8297 /* ------------------------------------------------------------------ */
8301 /* ------------------------------------------------------------------ */
8308 b0-=8; /* -> true start of storage */
8309 j=(uInt *)b0; /* -> first four bytes */
8313 b-b0-8, (Int)b0);
8316 b-b0-8, (Int)b0, n);
8318 decAllocBytes-=n; /* account for storage */
8319 /* printf(" free -- dAB: %d (%d)\n", decAllocBytes, -n); */