Skip to content

Commit 44a2bb6

Browse files
tompngtompng
authored
Round result of sqrt and BigMath methods (#427)
* Round result of sqrt and other BigMath methods BigMath methods no longer returns extra digits which may be inaccurate * Remove workaround for removing extra digits from sqrt result in log calculation
1 parent 2066c20 commit 44a2bb6

File tree

5 files changed

+69
-75
lines changed

5 files changed

+69
-75
lines changed

‎lib/bigdecimal.rb‎

Lines changed: 14 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -195,12 +195,9 @@ def sqrt(prec)
195195
raise FloatDomainError, "sqrt of 'NaN'(Not a Number)" if nan?
196196
return self if zero?
197197

198-
limit = BigDecimal.limit.nonzero? if prec == 0
199-
200-
# BigDecimal#sqrt calculates at least n_significant_digits precision.
201-
# This feature maybe problematic for some cases.
202-
n_digits = n_significant_digits
203-
prec = [prec, n_digits].max
198+
if prec == 0
199+
prec = BigDecimal.limit.nonzero? || n_significant_digits + BigDecimal.double_fig
200+
end
204201

205202
ex = exponent / 2
206203
x = _decimal_shift(-2 * ex)
@@ -210,8 +207,7 @@ def sqrt(prec)
210207
precs.reverse_each do |p|
211208
y = y.add(x.div(y, p), p).div(2, p)
212209
end
213-
y = y.mult(1, limit) if limit
214-
y._decimal_shift(ex)
210+
y._decimal_shift(ex).mult(1, prec)
215211
end
216212
end
217213

@@ -241,46 +237,43 @@ def self.log(x, prec)
241237
return BigDecimal::Internal.infinity_computation_result if x.infinite?
242238
return BigDecimal(0) if x == 1
243239

240+
prec2 = prec + BigDecimal.double_fig
244241
BigDecimal.save_limit do
245242
BigDecimal.limit(0)
246243
if x > 10 || x < 0.1
247-
log10 = log(BigDecimal(10), prec)
244+
log10 = log(BigDecimal(10), prec2)
248245
exponent = x.exponent
249246
x = x._decimal_shift(-exponent)
250247
if x < 0.3
251248
x *= 10
252249
exponent -= 1
253250
end
254-
return log10 * exponent + log(x, prec)
251+
return (log10 * exponent).add(log(x, prec2), prec)
255252
end
256253

257254
x_minus_one_exponent = (x - 1).exponent
258-
prec += BigDecimal.double_fig
259255

260256
# log(x) = log(sqrt(sqrt(sqrt(sqrt(x))))) * 2**sqrt_steps
261-
sqrt_steps = [Integer.sqrt(prec) + 3 * x_minus_one_exponent, 0].max
257+
sqrt_steps = [Integer.sqrt(prec2) + 3 * x_minus_one_exponent, 0].max
262258

263259
lg2 = 0.3010299956639812
264-
prec2 = prec + [-x_minus_one_exponent, 0].max + (sqrt_steps * lg2).ceil
260+
sqrt_prec = prec2 + [-x_minus_one_exponent, 0].max + (sqrt_steps * lg2).ceil
265261

266262
sqrt_steps.times do
267-
x = x.sqrt(prec2)
268-
269-
# Workaround for https://github.com/ruby/bigdecimal/issues/354
270-
x = x.mult(1, prec2 + BigDecimal.double_fig)
263+
x = x.sqrt(sqrt_prec)
271264
end
272265

273266
# Taylor series for log(x) around 1
274267
# log(x) = -log((1 + X) / (1 - X)) where X = (x - 1) / (x + 1)
275268
# log(x) = 2 * (X + X**3 / 3 + X**5 / 5 + X**7 / 7 + ...)
276-
x = (x - 1).div(x + 1, prec2)
269+
x = (x - 1).div(x + 1, sqrt_prec)
277270
y = x
278-
x2 = x.mult(x, prec)
271+
x2 = x.mult(x, prec2)
279272
1.step do |i|
280-
n = prec + x.exponent - y.exponent + x2.exponent
273+
n = prec2 + x.exponent - y.exponent + x2.exponent
281274
break if n <= 0 || x.zero?
282275
x = x.mult(x2.round(n - x2.exponent), n)
283-
y = y.add(x.div(2 * i + 1, n), prec)
276+
y = y.add(x.div(2 * i + 1, n), prec2)
284277
end
285278

286279
y.mult(2 ** (sqrt_steps + 1), prec)

‎lib/bigdecimal/math.rb‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -106,8 +106,8 @@ def sin(x, prec)
106106
d = x1.div(z,m)
107107
y += d
108108
end
109-
y *= sign
110-
y < -1 ? BigDecimal("-1") : y > 1 ? BigDecimal("1") : y
109+
y = BigDecimal("1") if y > 1
110+
y.mult(sign, prec)
111111
end
112112

113113
# call-seq:
@@ -145,7 +145,7 @@ def cos(x, prec)
145145
#
146146
def tan(x, prec)
147147
BigDecimal::Internal.validate_prec(prec, :tan)
148-
sin(x, prec) / cos(x, prec)
148+
sin(x, prec + BigDecimal.double_fig).div(cos(x, prec + BigDecimal.double_fig), prec)
149149
end
150150

151151
# call-seq:
@@ -163,11 +163,11 @@ def atan(x, prec)
163163
BigDecimal::Internal.validate_prec(prec, :atan)
164164
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :atan)
165165
return BigDecimal::Internal.nan_computation_result if x.nan?
166-
pi = PI(prec)
166+
n = prec + BigDecimal.double_fig
167+
pi = PI(n)
167168
x = -x if neg = x < 0
168169
return pi.div(neg ? -2 : 2, prec) if x.infinite?
169-
return pi / (neg ? -4 : 4) if x.round(prec) == 1
170-
n = prec + BigDecimal.double_fig
170+
return pi.div(neg ? -4 : 4, prec) if x.round(prec) == 1
171171
x = BigDecimal("1").div(x, n) if inv = x > 1
172172
x = (-1 + sqrt(1 + x.mult(x, n), n)).div(x, n) if dbl = x > 0.5
173173
y = x
@@ -185,7 +185,7 @@ def atan(x, prec)
185185
y *= 2 if dbl
186186
y = pi / 2 - y if inv
187187
y = -y if neg
188-
y
188+
y.mult(1, prec)
189189
end
190190

191191
# call-seq:
@@ -230,7 +230,7 @@ def PI(prec)
230230
pi = pi + d
231231
k = k+two
232232
end
233-
pi
233+
pi.mult(1, prec)
234234
end
235235

236236
# call-seq:

‎test/bigdecimal/helper.rb‎

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -37,16 +37,23 @@ def under_gc_stress
3737
GC.stress = stress
3838
end
3939

40-
# Asserts that the calculation of the given block converges to some value
41-
# with precision specified by block parameter.
40+
# Asserts that +actual+ is calculated with exactly the given +precision+.
41+
# No extra digits are allowed. Only the last digit may differ at most by one.
42+
def assert_in_exact_precision(expected, actual, precision)
43+
expected = BigDecimal(expected)
44+
delta = BigDecimal(1)._decimal_shift(expected.exponent - precision)
45+
assert actual.n_significant_digits <= precision, "Too many significant digits: #{actual.n_significant_digits} > #{precision}"
46+
assert_in_delta(expected.mult(1, precision), actual, delta)
47+
end
4248

49+
# Asserts that the calculation of the given block converges to some value
50+
# with exactly the given +precision+.
4351
def assert_converge_in_precision(&block)
4452
expected = yield(200)
4553
[50, 100, 150].each do |n|
4654
value = yield(n)
47-
precision = 1 - (value.div(expected, expected.precision) - 1).exponent
4855
assert(value != expected, "Unable to estimate precision for exact value")
49-
assert(precision >= n, "Precision is not enough: #{precision} < #{n}")
56+
assert_in_exact_precision(expected, value, n)
5057
end
5158
end
5259
end

‎test/bigdecimal/test_bigdecimal.rb‎

Lines changed: 12 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -1439,7 +1439,7 @@ def test_sqrt_bigdecimal
14391439
x = BigDecimal("-" + (2**100).to_s)
14401440
assert_raise_with_message(FloatDomainError, "sqrt of negative value") { x.sqrt(1) }
14411441
x = BigDecimal((2**200).to_s)
1442-
assert_equal(2**100, x.sqrt(1))
1442+
assert_equal(BigDecimal(2 ** 100).mult(1, 1), x.sqrt(1))
14431443

14441444
assert_in_delta(BigDecimal("4.0000000000000000000125"), BigDecimal("16.0000000000000000001").sqrt(100), BigDecimal("1e-40"))
14451445

@@ -1457,37 +1457,22 @@ def test_sqrt_bigdecimal
14571457

14581458
sqrt2_300 = BigDecimal(2).sqrt(300)
14591459
(250..270).each do |prec|
1460-
sqrt_prec = prec + BigDecimal.double_fig - 1
1461-
assert_in_delta(sqrt2_300, BigDecimal(2).sqrt(prec), BigDecimal("1e#{-sqrt_prec}"))
1460+
assert_in_exact_precision(sqrt2_300, BigDecimal(2).sqrt(prec), prec)
14621461
end
14631462
end
14641463

14651464
def test_sqrt_5266
14661465
x = BigDecimal('2' + '0'*100)
1467-
assert_equal('0.14142135623730950488016887242096980785696718753769480731',
1468-
x.sqrt(56).to_s(56).split(' ')[0])
1469-
assert_equal('0.1414213562373095048801688724209698078569671875376948073',
1470-
x.sqrt(55).to_s(55).split(' ')[0])
1466+
assert_equal('0.14142135623730950488016887242096980785696718753769480732e51',
1467+
x.sqrt(56).to_s)
1468+
assert_equal('0.1414213562373095048801688724209698078569671875376948073e51',
1469+
x.sqrt(55).to_s)
14711470

14721471
x = BigDecimal('2' + '0'*200)
1473-
assert_equal('0.14142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462',
1474-
x.sqrt(110).to_s(110).split(' ')[0])
1475-
assert_equal('0.1414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846',
1476-
x.sqrt(109).to_s(109).split(' ')[0])
1477-
end
1478-
1479-
def test_sqrt_minimum_precision
1480-
x = BigDecimal((2**200).to_s)
1481-
assert_equal(2**100, x.sqrt(1))
1482-
1483-
x = BigDecimal('1' * 60 + '.' + '1' * 40)
1484-
assert_in_delta(BigDecimal('3' * 30 + '.' + '3' * 70), x.sqrt(1), BigDecimal('1e-70'))
1485-
1486-
x = BigDecimal('1' * 40 + '.' + '1' * 60)
1487-
assert_in_delta(BigDecimal('3' * 20 + '.' + '3' * 80), x.sqrt(1), BigDecimal('1e-80'))
1488-
1489-
x = BigDecimal('0.' + '0' * 50 + '1' * 100)
1490-
assert_in_delta(BigDecimal('0.' + '0' * 25 + '3' * 100), x.sqrt(1), BigDecimal('1e-125'))
1472+
assert_equal('0.14142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462e101',
1473+
x.sqrt(110).to_s)
1474+
assert_equal('0.1414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846e101',
1475+
x.sqrt(109).to_s)
14911476
end
14921477

14931478
def test_internal_use_decimal_shift
@@ -2430,9 +2415,8 @@ def test_BigMath_log_with_square_of_E
24302415

24312416
def test_BigMath_log_with_high_precision_case
24322417
e = BigDecimal('2.71828182845904523536028747135266249775724709369996')
2433-
e_3 = e.mult(e, 50).mult(e, 50)
2434-
log_3 = BigMath.log(e_3, 50)
2435-
assert_in_delta(3, log_3, 0.0000000000_0000000000_0000000000_0000000000_0000000001)
2418+
e_3 = e.mult(e, 60).mult(e, 60)
2419+
assert_in_exact_precision(3, BigMath.log(e_3, 50), 50)
24362420
end
24372421

24382422
def test_BigMath_log_with_42

‎test/bigdecimal/test_bigmath.rb‎

Lines changed: 24 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ class TestBigMath < Test::Unit::TestCase
1616
def test_pi
1717
assert_equal(
1818
BigDecimal("3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068"),
19-
PI(100).round(100)
19+
PI(100)
2020
)
2121
assert_converge_in_precision {|n| PI(n) }
2222
end
@@ -37,8 +37,8 @@ def test_sqrt
3737
assert_raise(FloatDomainError) {sqrt(BigDecimal("-1.0"), N)}
3838
assert_raise(FloatDomainError) {sqrt(NAN, N)}
3939
assert_raise(FloatDomainError) {sqrt(PINF, N)}
40-
assert_in_delta(SQRT2, sqrt(BigDecimal("2"), 100), BigDecimal("1e-100"))
41-
assert_in_delta(SQRT3, sqrt(BigDecimal("3"), 100), BigDecimal("1e-100"))
40+
assert_in_exact_precision(SQRT2, sqrt(BigDecimal("2"), 100), 100)
41+
assert_in_exact_precision(SQRT3, sqrt(BigDecimal("3"), 100), 100)
4242
assert_converge_in_precision {|n| sqrt(BigDecimal("2"), n) }
4343
assert_converge_in_precision {|n| sqrt(BigDecimal("2e-50"), n) }
4444
assert_converge_in_precision {|n| sqrt(BigDecimal("2e50"), n) }
@@ -56,9 +56,9 @@ def test_sin
5656
assert_in_delta(0.0, sin(PI(N) * 21, N))
5757
assert_in_delta(0.0, sin(PI(N) * 30, N))
5858
assert_in_delta(-1.0, sin(PI(N) * BigDecimal("301.5"), N))
59-
assert_in_delta(BigDecimal('0.5'), sin(PI(100) / 6, 100), BigDecimal("1e-100"))
60-
assert_in_delta(SQRT3 / 2, sin(PI(100) / 3, 100), BigDecimal("1e-100"))
61-
assert_in_delta(SQRT2 / 2, sin(PI(100) / 4, 100), BigDecimal("1e-100"))
59+
assert_in_exact_precision(BigDecimal('0.5'), sin(PI(100) / 6, 100), 100)
60+
assert_in_exact_precision(SQRT3 / 2, sin(PI(100) / 3, 100), 100)
61+
assert_in_exact_precision(SQRT2 / 2, sin(PI(100) / 4, 100), 100)
6262
assert_converge_in_precision {|n| sin(BigDecimal("1"), n) }
6363
assert_converge_in_precision {|n| sin(BigDecimal("1e50"), n) }
6464
assert_converge_in_precision {|n| sin(BigDecimal("1e-30"), n) }
@@ -80,9 +80,9 @@ def test_cos
8080
assert_in_delta(-1.0, cos(PI(N) * 21, N))
8181
assert_in_delta(1.0, cos(PI(N) * 30, N))
8282
assert_in_delta(0.0, cos(PI(N) * BigDecimal("301.5"), N))
83-
assert_in_delta(BigDecimal('0.5'), cos(PI(100) / 3, 100), BigDecimal("1e-100"))
84-
assert_in_delta(SQRT3 / 2, cos(PI(100) / 6, 100), BigDecimal("1e-100"))
85-
assert_in_delta(SQRT2 / 2, cos(PI(100) / 4, 100), BigDecimal("1e-100"))
83+
assert_in_exact_precision(BigDecimal('0.5'), cos(PI(100) / 3, 100), 100)
84+
assert_in_exact_precision(SQRT3 / 2, cos(PI(100) / 6, 100), 100)
85+
assert_in_exact_precision(SQRT2 / 2, cos(PI(100) / 4, 100), 100)
8686
assert_converge_in_precision {|n| cos(BigDecimal("1"), n) }
8787
assert_converge_in_precision {|n| cos(BigDecimal("1e50"), n) }
8888
assert_converge_in_precision {|n| cos(BigDecimal(PI(50) / 2), n) }
@@ -100,14 +100,18 @@ def test_tan
100100
assert_in_delta(0.0, tan(-PI(N), N))
101101
assert_in_delta(-1.0, tan(-PI(N) / 4, N))
102102
assert_in_delta(-sqrt(BigDecimal(3), N), tan(-PI(N) / 3, N))
103+
assert_in_exact_precision(SQRT3, tan(PI(100) / 3, 100), 100)
104+
assert_converge_in_precision {|n| tan(1, n) }
105+
assert_converge_in_precision {|n| tan(BigMath::PI(50) / 2, n) }
106+
assert_converge_in_precision {|n| tan(BigMath::PI(50), n) }
103107
end
104108

105109
def test_atan
106110
assert_equal(0.0, atan(BigDecimal("0.0"), N))
107111
assert_in_delta(Math::PI/4, atan(BigDecimal("1.0"), N))
108112
assert_in_delta(Math::PI/6, atan(sqrt(BigDecimal("3.0"), N) / 3, N))
109113
assert_in_delta(Math::PI/2, atan(PINF, N))
110-
assert_in_delta(PI(100) / 3, atan(SQRT3, 100), BigDecimal("1e-100"))
114+
assert_in_exact_precision(PI(100) / 3, atan(SQRT3, 100), 100)
111115
assert_equal(BigDecimal("0.823840753418636291769355073102514088959345624027952954058347023122539489"),
112116
atan(BigDecimal("1.08"), 72).round(72), '[ruby-dev:41257]')
113117
assert_converge_in_precision {|n| atan(BigDecimal("2"), n)}
@@ -120,8 +124,11 @@ def test_exp
120124
assert_in_epsilon(Math.exp(x), BigMath.exp(BigDecimal(x, 0), N))
121125
end
122126
assert_equal(1, BigMath.exp(BigDecimal("0"), N))
123-
assert_in_epsilon(BigDecimal("4.48168907033806482260205546011927581900574986836966705677265008278593667446671377298105383138245339138861635065183019577"),
124-
BigMath.exp(BigDecimal("1.5"), 100), BigDecimal("1e-100"))
127+
assert_in_exact_precision(
128+
BigDecimal("4.48168907033806482260205546011927581900574986836966705677265008278593667446671377298105383138245339138861635065183019577"),
129+
BigMath.exp(BigDecimal("1.5"), 100),
130+
100
131+
)
125132
assert_converge_in_precision {|n| BigMath.exp(BigDecimal("1"), n) }
126133
assert_converge_in_precision {|n| BigMath.exp(BigDecimal("-2"), n) }
127134
assert_converge_in_precision {|n| BigMath.exp(BigDecimal("-34"), n) }
@@ -132,8 +139,11 @@ def test_exp
132139
def test_log
133140
assert_equal(0, BigMath.log(BigDecimal("1.0"), 10))
134141
assert_in_epsilon(Math.log(10)*1000, BigMath.log(BigDecimal("1e1000"), 10))
135-
assert_in_epsilon(BigDecimal("2.3025850929940456840179914546843642076011014886287729760333279009675726096773524802359972050895982983419677840422862"),
136-
BigMath.log(BigDecimal("10"), 100), BigDecimal("1e-100"))
142+
assert_in_exact_precision(
143+
BigDecimal("2.3025850929940456840179914546843642076011014886287729760333279009675726096773524802359972050895982983419677840422862"),
144+
BigMath.log(BigDecimal("10"), 100),
145+
100
146+
)
137147
assert_converge_in_precision {|n| BigMath.log(BigDecimal("2"), n) }
138148
assert_converge_in_precision {|n| BigMath.log(BigDecimal("1e-30") + 1, n) }
139149
assert_converge_in_precision {|n| BigMath.log(BigDecimal("1e-30"), n) }

0 commit comments

Comments
 (0)