ruby-changes:17291
From: mrkn <ko1@a...>
Date: Sun, 19 Sep 2010 05:23:56 +0900 (JST)
Subject: [ruby-changes:17291] Ruby:r29293 (trunk): * ext/bigdecimal/bigdecimal.c: fix rounding algorithms for half-down
mrkn 2010-09-19 05:19:38 +0900 (Sun, 19 Sep 2010) New Revision: 29293 http://svn.ruby-lang.org/cgi-bin/viewvc.cgi?view=rev&revision=29293 Log: * ext/bigdecimal/bigdecimal.c: fix rounding algorithms for half-down and half-even. This change is based on the patch created by Matthew Willson, the reporter of this bug. [Bug #3803] [ruby-core:32136] * test/bigdecimal/test_bigdecimal.rb: add tests for above changes. Modified files: trunk/ChangeLog trunk/ext/bigdecimal/bigdecimal.c trunk/test/bigdecimal/test_bigdecimal.rb Index: ChangeLog =================================================================== --- ChangeLog (revision 29292) +++ ChangeLog (revision 29293) @@ -1,3 +1,11 @@ +Sun Sep 19 05:14:35 2010 Kenta Murata <mrkn@m...> + + * ext/bigdecimal/bigdecimal.c: fix rounding algorithms for half-down + and half-even. This change is based on the patch created by Matthew + Willson, the reporter of this bug. [Bug #3803] [ruby-core:32136] + + * test/bigdecimal/test_bigdecimal.rb: add tests for above changes. + Sat Sep 18 20:09:51 2010 Tanaka Akira <akr@f...> * ext/pathname/pathname.c (path_each_entry): Pathname#each_entry Index: ext/bigdecimal/bigdecimal.c =================================================================== --- ext/bigdecimal/bigdecimal.c (revision 29292) +++ ext/bigdecimal/bigdecimal.c (revision 29293) @@ -308,9 +308,9 @@ * * ROUND_UP:: round away from zero * ROUND_DOWN:: round towards zero (truncate) - * ROUND_HALF_UP:: round up if the appropriate digit >= 5, otherwise truncate (default) - * ROUND_HALF_DOWN:: round up if the appropriate digit >= 6, otherwise truncate - * ROUND_HALF_EVEN:: round towards the even neighbor (Banker's rounding) + * ROUND_HALF_UP:: round towards the nearest neighbor, unless both neighbors are equidistant, in which case round away from zero. (default) + * ROUND_HALF_DOWN:: round towards the nearest neighbor, unless both neighbors are equidistant, in which case round towards zero. + * ROUND_HALF_EVEN:: round towards the nearest neighbor, unless both neighbors are equidistant, in which case round towards the even neighbor (Banker's rounding) * ROUND_CEILING:: round towards positive infinity (ceil) * ROUND_FLOOR:: round towards negative infinity (floor) * @@ -4559,8 +4559,9 @@ */ { /* fracf: any positive digit under rounding position? */ + /* fracf_1further: any positive digits under one further than the rounding position? */ /* exptoadd: number of digits needed to compensate negative nf */ - int fracf; + int fracf, fracf_1further; ssize_t n,i,ix,ioffset, exptoadd; BDIGIT v, shifter; BDIGIT div; @@ -4573,62 +4574,111 @@ VpSetZero(y,VpGetSign(y)); /* truncate everything */ return 0; } - exptoadd = -nf; - nf = 0; + exptoadd = -nf; + nf = 0; } - /* ix: x->fraq[ix] contains round position */ ix = nf / (ssize_t)BASE_FIG; if ((size_t)ix >= y->Prec) return 0; /* rounding position too right(small). */ - ioffset = nf - ix*(ssize_t)BASE_FIG; - v = y->frac[ix]; - /* drop digits after pointed digit */ + ioffset = nf - ix*(ssize_t)BASE_FIG; n = (ssize_t)BASE_FIG - ioffset - 1; for (shifter=1,i=0; i<n; ++i) shifter *= 10; + + /* so the representation used (in y->frac) is an array of BDIGIT, where + each BDIGIT contains a value between 0 and BASE-1, consisting of BASE_FIG + decimal places. + + (that numbers of decimal places are typed as ssize_t is somewhat confusing) + + nf is now position (in decimal places) of the digit from the start of + the array. + ix is the position (in BDIGITS) of the BDIGIT containing the decimal digit, + from the start of the array. + v is the value of this BDIGIT + ioffset is the number of extra decimal places along of this decimal digit + within v. + n is the number of decimal digits remaining within v after this decimal digit + shifter is 10**n, + v % shifter are the remaining digits within v + v % (shifter * 10) are the digit together with the remaining digits within v + v / shifter are the digit's predecessors together with the digit + div = v / shifter / 10 is just the digit's precessors + (v / shifter) - div*10 is just the digit, which is what v ends up being reassigned to. + */ + fracf = (v % (shifter * 10) > 0); + fracf_1further = ((v % shifter) > 0); + v /= shifter; div = v / 10; v = v - div*10; - if (fracf == 0) { - for (i=ix+1; (size_t)i < y->Prec; i++) { - if (y->frac[i] % BASE) { - fracf = 1; - break; - } - } + /* now v is just the digit required. + now fracf is whether the digit or any of the remaining digits within v are non-zero + now fracf_1further is whether any of the remaining digits within v are non-zero + */ + + /* now check all the remaining BDIGITS for zero-ness a whole BDIGIT at a time. + if we spot any non-zeroness, that means that we foudn a positive digit under + rounding position, and we also found a positive digit under one further than + the rounding position, so both searches (to see if any such non-zero digit exists) + can stop */ + + for (i=ix+1; (size_t)i < y->Prec; i++) { + if (y->frac[i] % BASE) { + fracf = fracf_1further = 1; + break; + } } + + /* now fracf = does any positive digit exist under the rounding position? + now fracf_1further = does any positive digit exist under one further than the + rounding position? + now v = the first digit under the rounding position */ + + /* drop digits after pointed digit */ memset(y->frac+ix+1, 0, (y->Prec - (ix+1)) * sizeof(BDIGIT)); + switch(f) { case VP_ROUND_DOWN: /* Truncate */ break; case VP_ROUND_UP: /* Roundup */ if (fracf) ++div; - break; - case VP_ROUND_HALF_UP: /* Round half up */ + break; + case VP_ROUND_HALF_UP: if (v>=5) ++div; break; - case VP_ROUND_HALF_DOWN: /* Round half down */ - if (v>=6) ++div; + case VP_ROUND_HALF_DOWN: + if (v > 5 || (v == 5 && fracf_1further)) ++div; break; - case VP_ROUND_CEIL: /* ceil */ + case VP_ROUND_CEIL: if (fracf && (VpGetSign(y)>0)) ++div; break; - case VP_ROUND_FLOOR: /* floor */ + case VP_ROUND_FLOOR: if (fracf && (VpGetSign(y)<0)) ++div; break; case VP_ROUND_HALF_EVEN: /* Banker's rounding */ - if (v>5) ++div; - else if (v==5) { - if ((size_t)i == (BASE_FIG-1)) { - if (ix && (y->frac[ix-1]%2)) ++div; - } + if (v > 5) ++div; + else if (v == 5) { + if (fracf_1further) { + ++div; + } else { - if (div%2) ++div; - } - } - break; + if (ioffset == 0) { + /* v is the first decimal digit of its BDIGIT; + need to grab the previous BDIGIT if present + to check for evenness of the previous decimal + digit (which is same as that of the BDIGIT since + base 10 has a factor of 2) */ + if (ix && (y->frac[ix-1] % 2)) ++div; + } + else { + if (div % 2) ++div; + } + } + } + break; } for (i=0; i<=n; ++i) div *= 10; if (div>=BASE) { @@ -4694,6 +4744,8 @@ return VpLeftRound(c, (int)VpGetRoundMode(), (ssize_t)ix); } +/* If I understand correctly, this is only ever used to round off the final decimal + digit of precision */ static void VpInternalRound(Real *c, size_t ixDigit, BDIGIT vPrev, BDIGIT v) { @@ -4701,36 +4753,41 @@ unsigned short const rounding_mode = VpGetRoundMode(); - if(VpLimitRound(c,ixDigit)) return; - if(!v) return; + if (VpLimitRound(c, ixDigit)) return; + if (!v) return; v /= BASE1; switch (rounding_mode) { case VP_ROUND_DOWN: - break; + break; case VP_ROUND_UP: - if(v) f = 1; - break; + if (v) f = 1; + break; case VP_ROUND_HALF_UP: - if(v >= 5) f = 1; - break; + if (v >= 5) f = 1; + break; case VP_ROUND_HALF_DOWN: - if(v >= 6) f = 1; - break; - case VP_ROUND_CEIL: /* ceil */ - if(v && (VpGetSign(c)>0)) f = 1; - break; - case VP_ROUND_FLOOR: /* floor */ - if(v && (VpGetSign(c)<0)) f = 1; - break; + /* this is ok - because this is the last digit of precision, + the case where v == 5 and some further digits are nonzero + will never occur */ + if (v >= 6) f = 1; + break; + case VP_ROUND_CEIL: + if (v && (VpGetSign(c) > 0)) f = 1; + break; + case VP_ROUND_FLOOR: + if (v && (VpGetSign(c) < 0)) f = 1; + break; case VP_ROUND_HALF_EVEN: /* Banker's rounding */ - if(v>5) f = 1; - else if(v==5 && vPrev%2) f = 1; - break; + /* as per VP_ROUND_HALF_DOWN, because this is the last digit of precision, + there is no case to worry about where v == 5 and some further digits are nonzero */ + if (v > 5) f = 1; + else if (v == 5 && vPrev % 2) f = 1; + break; } - if(f) { - VpRdup(c,ixDigit); /* round up */ - VpNmlz(c); + if (f) { + VpRdup(c, ixDigit); + VpNmlz(c); } } Index: test/bigdecimal/test_bigdecimal.rb =================================================================== --- test/bigdecimal/test_bigdecimal.rb (revision 29292) +++ test/bigdecimal/test_bigdecimal.rb (revision 29293) @@ -612,6 +612,18 @@ assert_equal(3, x.round(0, BigDecimal::ROUND_CEILING)) assert_equal(2, x.round(0, BigDecimal::ROUND_FLOOR)) assert_raise(TypeError) { x.round(0, 256) } + + 15.times do |n| + x = BigDecimal.new("5#{'0'*n}1") + assert_equal(10**(n+2), x.round(-(n+2), BigDecimal::ROUND_HALF_DOWN)) + assert_equal(10**(n+2), x.round(-(n+2), BigDecimal::ROUND_HALF_EVEN)) + x = BigDecimal.new("0.5#{'0'*n}1") + assert_equal(1, x.round(0, BigDecimal::ROUND_HALF_DOWN)) + assert_equal(1, x.round(0, BigDecimal::ROUND_HALF_EVEN)) + x = BigDecimal.new("-0.5#{'0'*n}1") + assert_equal(-1, x.round(0, BigDecimal::ROUND_HALF_DOWN)) + assert_equal(-1, x.round(0, BigDecimal::ROUND_HALF_EVEN)) + end end def test_truncate -- ML: ruby-changes@q... Info: http://www.atdot.net/~ko1/quickml/