Stefan Krah's "beta quality" fixes for issues that have been discovered in arbitrary-length (decNumber) arithmetic when one or both of the input operands has a precision that is larger than the context precision. Two of these fix earlier changesets, so they have to be applied on top of each other. ----- # HG changeset patch # User Stefan Krah # Date 1391612403 -3600 # Wed Feb 05 16:00:03 2014 +0100 # Node ID 9251a4072b87bcdb70f859171f3d5c8629e70617 # Parent 2d30038beaa4f28709730028c38bc392cd43f8aa Fix rounding and flags for add/subtract/fma with overlong inputs. In decAddOp(), the condition "rhs->digits+padding > lhs->digits+reqdigits+1", which triggers code to avoid huge shifts, is not restrictive enough for overlong inputs. Consider for example 100000005e6 - 5000000, with a context precision of 4 and ROUND_CEILING: | rhs | 100000005000000 5000000 | lhs | expected: 1.000E+14 Rounded calculated: 1.001E+14 Inexact Rounded diff --git a/decNumber.c b/decNumber.c --- a/decNumber.c +++ b/decNumber.c @@ -3805,6 +3805,7 @@ decNumber *alloclhs=NULL; // non-NULL if rounded lhs allocated decNumber *allocrhs=NULL; // .., rhs #endif + decNumber dtiny; // lhs, adjusted to avoid a huge shift Int rhsshift; // working shift (in Units) Int maxdigits; // longest logical length Int mult; // multiplier @@ -3969,12 +3970,7 @@ } // Now align (pad) the lhs or rhs so they can be added or - // subtracted, as necessary. If one number is much larger than - // the other (that is, if in plain form there is a least one - // digit between the lowest digit of one and the highest of the - // other) padding with up to DIGITS-1 trailing zeros may be - // needed; then apply rounding (as exotic rounding modes may be - // affected by the residue). + // subtracted, as necessary. rhsshift=0; // rhs shift to left (padding) in Units bits=lhs->bits; // assume sign is that of LHS mult=1; // likely multiplier @@ -3984,35 +3980,26 @@ // some padding needed; always pad the RHS, as any required // padding can then be effected by a simple combination of // shifts and a multiply - Flag swapped=0; + Int exponent; // new exponent of LHS (if adjusted) if (padding<0) { // LHS needs the padding const decNumber *t; padding=-padding; // will be +ve bits=(uByte)(rhs->bits^negate); // assumed sign is now that of RHS t=lhs; lhs=rhs; rhs=t; - swapped=1; } - // If, after pad, rhs would be longer than lhs by digits+1 or - // more then lhs cannot affect the answer, except as a residue, - // so only need to pad up to a length of DIGITS+1. - if (rhs->digits+padding > lhs->digits+reqdigits+1) { - // The RHS is sufficient - // for residue use the relative sign indication... - Int shift=reqdigits-rhs->digits; // left shift needed - residue=1; // residue for rounding - if (diffsign) residue=-residue; // signs differ - // copy, shortening if necessary - decCopyFit(res, rhs, set, &residue, status); - // if it was already shorter, then need to pad with zeros - if (shift>0) { - res->digits=decShiftToMost(res->lsu, res->digits, shift); - res->exponent-=shift; // adjust the exponent. - } - // flip the result sign if unswapped and rhs was negated - if (!swapped) res->bits^=negate; - decFinish(res, set, &residue, status); // done - break;} + exponent = rhs->exponent-1; + exponent += (rhs->digits>reqdigits) ? 0 : rhs->digits-reqdigits-1; + if (lhs->digits+lhs->exponent-1 < exponent) { + // Adjust lhs and padding to avoid huge shifts. + dtiny.bits=lhs->bits; + dtiny.exponent=exponent; + dtiny.digits=1; + dtiny.lsu[0]=1; + lhs=&dtiny; + padding=rhs->exponent-exponent; + // Fall through to add/subtract the modified lhs. + } // LHS digits may affect result rhsshift=D2U(padding+1)-1; // this much by Unit shift .. # HG changeset patch # User Stefan Krah # Date 1391628142 -3600 # Wed Feb 05 20:22:22 2014 +0100 # Node ID cc70d6dff9e9e5c4b6aea8d6bfc5ec5dd31320ce # Parent 9251a4072b87bcdb70f859171f3d5c8629e70617 Fix incorrect results with overlong input in NextMinus/NextPlus. diff --git a/decNumber.c b/decNumber.c --- a/decNumber.c +++ b/decNumber.c @@ -1668,10 +1668,18 @@ // there is no status to set return res; } - decNumberZero(&dtiny); // start with 0 - dtiny.lsu[0]=1; // make number that is .. - dtiny.exponent=DEC_MIN_EMIN-1; // .. smaller than tiniest + // Apply the context. If the result is inexact, we are done. workset.round=DEC_ROUND_FLOOR; + workset.status=0; + decNumberPlus(res, rhs, &workset); + if (workset.status&(DEC_Inexact|DEC_NaNs)) { + set->status |= (workset.status&DEC_NaNs); + return res; + } + // Applying the context was exact. Subtract tiny quantity and round. + decNumberZero(&dtiny); // start with 0 + dtiny.lsu[0]=1; + dtiny.exponent=set->emin-set->digits; decAddOp(res, rhs, &dtiny, &workset, DECNEG, &status); status&=DEC_Invalid_operation|DEC_sNaN; // only sNaN Invalid please if (status!=0) decStatus(res, status, set); @@ -1705,10 +1713,18 @@ // there is no status to set return res; } - decNumberZero(&dtiny); // start with 0 - dtiny.lsu[0]=1; // make number that is .. - dtiny.exponent=DEC_MIN_EMIN-1; // .. smaller than tiniest + // Apply the context. If the result is inexact, we are done. workset.round=DEC_ROUND_CEILING; + workset.status=0; + decNumberPlus(res, rhs, &workset); + if (workset.status&(DEC_Inexact|DEC_NaNs)) { + set->status |= (workset.status|DEC_NaNs); + return res; + } + // Applying the context was exact. Add tiny quantity and round. + decNumberZero(&dtiny); // start with 0 + dtiny.lsu[0]=1; + dtiny.exponent=set->emin-set->digits; decAddOp(res, rhs, &dtiny, &workset, 0, &status); status&=DEC_Invalid_operation|DEC_sNaN; // only sNaN Invalid please if (status!=0) decStatus(res, status, set); # HG changeset patch # User Stefan Krah # Date 1391628552 -3600 # Wed Feb 05 20:29:12 2014 +0100 # Node ID 1a08ae1cf0c88f76bc3ab7ea6023c8a0ea7aab7f # Parent cc70d6dff9e9e5c4b6aea8d6bfc5ec5dd31320ce Fix incorrect power result with overlong input. This is a preliminary fix that should be audited in depth (though it looks "probably correct"). diff --git a/decNumber.c b/decNumber.c --- a/decNumber.c +++ b/decNumber.c @@ -2132,6 +2132,9 @@ aset.round=DEC_ROUND_HALF_EVEN; // internally use balanced // calculate the working DIGITS aset.digits=reqdigits+(rhs->digits+rhs->exponent)+2; + aset.emax=DEC_MAX_MATH; // usual bounds + aset.emin=-DEC_MAX_MATH; // .. + aset.clamp=0; // and no concrete format #if DECSUBSET if (!set->extended) aset.digits--; // use classic precision #endif # HG changeset patch # User Stefan Krah # Date 1391629669 -3600 # Wed Feb 05 20:47:49 2014 +0100 # Node ID bc7c8815d2460e5cf47381c7d85b448a64d4e0d4 # Parent 1a08ae1cf0c88f76bc3ab7ea6023c8a0ea7aab7f Conservative fix for an infinite loop in ln with overlong input. The fix has not been formally analyzed. diff --git a/decNumber.c b/decNumber.c --- a/decNumber.c +++ b/decNumber.c @@ -5729,6 +5729,10 @@ aset.emax=set->emax; aset.emin=set->emin; aset.clamp=0; // no concrete format + if (rhs->digits > set->digits) { + aset.emax=DEC_MAX_MATH*2; + aset.emin=-DEC_MAX_MATH*2; + } // set up a context to be used for the multiply and subtract bset=aset; bset.emax=DEC_MAX_MATH*2; // use double bounds for the # HG changeset patch # User Stefan Krah # Date 1391629728 -3600 # Wed Feb 05 20:48:48 2014 +0100 # Node ID c6f22597a4fb4038511d321950b2e66cb78158da # Parent bc7c8815d2460e5cf47381c7d85b448a64d4e0d4 Fix typo in the NextPlus patch. diff --git a/decNumber.c b/decNumber.c --- a/decNumber.c +++ b/decNumber.c @@ -1718,7 +1718,7 @@ workset.status=0; decNumberPlus(res, rhs, &workset); if (workset.status&(DEC_Inexact|DEC_NaNs)) { - set->status |= (workset.status|DEC_NaNs); + set->status |= (workset.status&DEC_NaNs); return res; } // Applying the context was exact. Add tiny quantity and round. # HG changeset patch # User Stefan Krah # Date 1391630917 -3600 # Wed Feb 05 21:08:37 2014 +0100 # Node ID cda06c664d57ebc0173a2c7867adcf9e595b3934 # Parent c6f22597a4fb4038511d321950b2e66cb78158da Use DEC_MAX_EMAX / DEC_MIN_EMIN, since this is the integer exponent path. diff --git a/decNumber.c b/decNumber.c --- a/decNumber.c +++ b/decNumber.c @@ -2132,8 +2132,8 @@ aset.round=DEC_ROUND_HALF_EVEN; // internally use balanced // calculate the working DIGITS aset.digits=reqdigits+(rhs->digits+rhs->exponent)+2; - aset.emax=DEC_MAX_MATH; // usual bounds - aset.emin=-DEC_MAX_MATH; // .. + aset.emax=DEC_MAX_EMAX; // usual bounds + aset.emin=DEC_MIN_EMIN; // .. aset.clamp=0; // and no concrete format #if DECSUBSET if (!set->extended) aset.digits--; // use classic precision # HG changeset patch # User Stefan Krah # Date 1391631931 -3600 # Wed Feb 05 21:25:31 2014 +0100 # Node ID 0803e97f53e1769e5db5ae1e94ad83c168502e88 # Parent cda06c664d57ebc0173a2c7867adcf9e595b3934 Decapitate overlong input before use in rotate. diff --git a/decNumber.c b/decNumber.c --- a/decNumber.c +++ b/decNumber.c @@ -2489,6 +2489,7 @@ status=DEC_Invalid_operation; else { // rhs is OK decNumberCopy(res, lhs); + if (res->digits>set->digits) decDecap(res, res->digits-set->digits); // convert -ve rotate to equivalent positive rotation if (rotate<0) rotate=set->digits+rotate; if (rotate!=0 && rotate!=set->digits // zero or full rotation # HG changeset patch # User Stefan Krah # Date 1391687067 -3600 # Thu Feb 06 12:44:27 2014 +0100 # Node ID b154a8783fa6b234fd83112b411d2f270d4dfdb8 # Parent 0803e97f53e1769e5db5ae1e94ad83c168502e88 Change scaleb to use the context precision for overlong inputs. diff --git a/decNumber.c b/decNumber.c --- a/decNumber.c +++ b/decNumber.c @@ -2660,6 +2660,7 @@ else res->exponent=DEC_MAX_EMAX+1; } residue=0; + decCopyFit(res, res, set, &residue, &status); decFinalize(res, set, &residue, &status); // final check } // finite LHS } // rhs OK # HG changeset patch # User Stefan Krah # Date 1391690745 -3600 # Thu Feb 06 13:45:45 2014 +0100 # Node ID 7b6724db1a1c2e777425ee5bd250cc9fd2dd49c8 # Parent b154a8783fa6b234fd83112b411d2f270d4dfdb8 log10: Handle payloads that exceed the context precision. diff --git a/decNumber.c b/decNumber.c --- a/decNumber.c +++ b/decNumber.c @@ -1408,6 +1408,12 @@ if (decCheckOperands(res, DECUNUSED, rhs, set)) return res; #endif + // Handle payloads that exceed the context precision. + if (rhs->bits&(DECNAN|DECSNAN)) { + decNumberPlus(res, rhs, set); + return res; + } + // Check restrictions; this is a math function; if not violated // then carry out the operation. if (!decCheckMath(rhs, set, &status)) do { // protect malloc # HG changeset patch # User Stefan Krah # Date 1391695572 -3600 # Thu Feb 06 15:06:12 2014 +0100 # Node ID 614907e5f165977ac45b03332f01bd3ba8941d80 # Parent 7b6724db1a1c2e777425ee5bd250cc9fd2dd49c8 decDivideOp: Always check the exponent if lhs==0 to ensure proper clamping. diff --git a/decNumber.c b/decNumber.c --- a/decNumber.c +++ b/decNumber.c @@ -4356,12 +4356,10 @@ else { #endif if (op&DIVIDE) { - residue=0; exponent=lhs->exponent-rhs->exponent; // ideal exponent decNumberCopy(res, lhs); // [zeros always fit] res->bits=bits; // sign as computed res->exponent=exponent; // exponent, too - decFinalize(res, set, &residue, status); // check exponent } else if (op&DIVIDEINT) { decNumberZero(res); // integer 0 @@ -4372,6 +4370,8 @@ decNumberCopy(res, lhs); // [zeros always fit] if (exponentexponent) res->exponent=exponent; // use lower } + residue=0; + decFinalize(res, set, &residue, status); // check exponent #if DECSUBSET } #endif