Skip to content

Rewrite BigDecimal#sqrt in ruby using Integer.sqrt #323

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
tompng wants to merge 6 commits into from
Closed

Rewrite BigDecimal#sqrt in ruby using Integer.sqrt #323

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b7dfb3b
Rewrite BigDecimal#sqrt in ruby using Integer.sqrt
tompng May 9, 2025
b1826f2
Add a fast path for sqrt when sqrt is an exact value
tompng May 18, 2025
03400c1
Use BigDecimal("1e#{n}") instead of ten**n to make calculation faster
tompng May 19, 2025
7338d17
Tweak fast path exponent calculation and add more tests
tompng May 20, 2025
7de5f52
Remove unused VpConstPt5 used in sqrt before
tompng Jun 17, 2025
2ee3a3b
pend sqrt test in older version of truffleruby
tompng Jun 18, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
202 changes: 0 additions & 202 deletions ext/bigdecimal/bigdecimal.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2355,31 +2355,6 @@ BigDecimal_abs(VALUE self)
return CheckGetValue(c);
}

/* call-seq:
* sqrt(n)
*
* Returns the square root of the value.
*
* Result has at least n significant digits.
*/
static VALUE
BigDecimal_sqrt(VALUE self, VALUE nFig)
{
ENTER(5);
BDVALUE c, a;
size_t mx, n;

GUARD_OBJ(a, GetBDValueMust(self));
mx = a.real->Prec * (VpBaseFig() + 1);

n = check_int_precision(nFig);
n += VpDblFig() + VpBaseFig();
if (mx <= n) mx = n;
GUARD_OBJ(c, NewZeroWrapLimited(1, mx));
VpSqrt(c.real, a.real);
return CheckGetValue(c);
}

/* Return the integer part of the number, as a BigDecimal.
*/
static VALUE
Expand Down Expand Up @@ -4577,7 +4552,6 @@ Init_bigdecimal(void)
rb_define_method(rb_cBigDecimal, "dup", BigDecimal_clone, 0);
rb_define_method(rb_cBigDecimal, "to_f", BigDecimal_to_f, 0);
rb_define_method(rb_cBigDecimal, "abs", BigDecimal_abs, 0);
rb_define_method(rb_cBigDecimal, "sqrt", BigDecimal_sqrt, 1);
rb_define_method(rb_cBigDecimal, "fix", BigDecimal_fix, 0);
rb_define_method(rb_cBigDecimal, "round", BigDecimal_round, -1);
rb_define_method(rb_cBigDecimal, "frac", BigDecimal_frac, 0);
Expand Down Expand Up @@ -4655,9 +4629,6 @@ static int gfDebug = 1; /* Debug switch */
#endif /* BIGDECIMAL_DEBUG */

static Real *VpConstOne; /* constant 1.0 */
static Real *VpConstPt5; /* constant 0.5 */
#define maxnr 100UL /* Maximum iterations for calculating sqrt. */
/* used in VpSqrt() */

enum op_sw {
OP_SW_ADD = 1, /* + */
Expand Down Expand Up @@ -5063,11 +5034,6 @@ VpInit(DECDIG BaseVal)
/* Const 1.0 */
VpConstOne = NewOneNolimit(1, 1);

/* Const 0.5 */
VpConstPt5 = NewOneNolimit(1, 1);
VpConstPt5->exponent = 0;
VpConstPt5->frac[0] = 5*BASE1;

#ifdef BIGDECIMAL_DEBUG
gnAlloc = 0;
#endif /* BIGDECIMAL_DEBUG */
Expand Down Expand Up @@ -6892,174 +6858,6 @@ VpVtoD(double *d, SIGNED_VALUE *e, Real *m)
return f;
}

/*
* m <- d
*/
VP_EXPORT void
VpDtoV(Real *m, double d)
{
size_t ind_m, mm;
SIGNED_VALUE ne;
DECDIG i;
double val, val2;

if (isnan(d)) {
VpSetNaN(m);
goto Exit;
}
if (isinf(d)) {
if (d > 0.0) VpSetPosInf(m);
else VpSetNegInf(m);
goto Exit;
}

if (d == 0.0) {
VpSetZero(m, 1);
goto Exit;
}
val = (d > 0.) ? d : -d;
ne = 0;
if (val >= 1.0) {
while (val >= 1.0) {
val /= (double)BASE;
++ne;
}
}
else {
val2 = 1.0 / (double)BASE;
while (val < val2) {
val *= (double)BASE;
--ne;
}
}
/* Now val = 0.xxxxx*BASE**ne */

mm = m->MaxPrec;
memset(m->frac, 0, mm * sizeof(DECDIG));
for (ind_m = 0; val > 0.0 && ind_m < mm; ind_m++) {
val *= (double)BASE;
i = (DECDIG)val;
val -= (double)i;
m->frac[ind_m] = i;
}
if (ind_m >= mm) ind_m = mm - 1;
VpSetSign(m, (d > 0.0) ? 1 : -1);
m->Prec = ind_m + 1;
m->exponent = ne;

VpInternalRound(m, 0, (m->Prec > 0) ? m->frac[m->Prec-1] : 0,
(DECDIG)(val*(double)BASE));

Exit:
return;
}

/*
* y = SQRT(x), y*y - x =>0
*/
VP_EXPORT int
VpSqrt(Real *y, Real *x)
{
Real *f = NULL;
Real *r = NULL;
size_t y_prec;
SIGNED_VALUE n, e;
ssize_t nr;
double val;

/* Zero or +Infinity ? */
if (VpIsZero(x) || VpIsPosInf(x)) {
VpAsgn(y,x,1);
goto Exit;
}

/* Negative ? */
if (BIGDECIMAL_NEGATIVE_P(x)) {
VpSetNaN(y);
return VpException(VP_EXCEPTION_OP, "sqrt of negative value", 0);
}

/* NaN ? */
if (VpIsNaN(x)) {
VpSetNaN(y);
return VpException(VP_EXCEPTION_OP, "sqrt of 'NaN'(Not a Number)", 0);
}

/* One ? */
if (VpIsOne(x)) {
VpSetOne(y);
goto Exit;
}

n = (SIGNED_VALUE)y->MaxPrec;
if (x->MaxPrec > (size_t)n) n = (ssize_t)x->MaxPrec;

/* allocate temporally variables */
/* TODO: reconsider MaxPrec of f and r */
f = NewOneNolimit(1, y->MaxPrec * (BASE_FIG + 2));
r = NewOneNolimit(1, (n + n) * (BASE_FIG + 2));

nr = 0;
y_prec = y->MaxPrec;

VpVtoD(&val, &e, x); /* val <- x */
e /= (SIGNED_VALUE)BASE_FIG;
n = e / 2;
if (e - n * 2 != 0) {
val /= BASE;
n = (e + 1) / 2;
}
VpDtoV(y, sqrt(val)); /* y <- sqrt(val) */
y->exponent += n;
n = (SIGNED_VALUE)roomof(BIGDECIMAL_DOUBLE_FIGURES, BASE_FIG);
y->MaxPrec = Min((size_t)n , y_prec);
f->MaxPrec = y->MaxPrec + 1;
n = (SIGNED_VALUE)(y_prec * BASE_FIG);
if (n > (SIGNED_VALUE)maxnr) n = (SIGNED_VALUE)maxnr;

/*
* Perform: y_{n+1} = (y_n - x/y_n) / 2
*/
do {
y->MaxPrec *= 2;
if (y->MaxPrec > y_prec) y->MaxPrec = y_prec;
f->MaxPrec = y->MaxPrec;
VpDivd(f, r, x, y); /* f = x/y */
VpAddSub(r, f, y, -1); /* r = f - y */
VpMult(f, VpConstPt5, r); /* f = 0.5*r */
if (y_prec == y->MaxPrec && VpIsZero(f))
goto converge;
VpAddSub(r, f, y, 1); /* r = y + f */
VpAsgn(y, r, 1); /* y = r */
} while (++nr < n);

#ifdef BIGDECIMAL_DEBUG
if (gfDebug) {
printf("ERROR(VpSqrt): did not converge within %ld iterations.\n", nr);
}
#endif /* BIGDECIMAL_DEBUG */
y->MaxPrec = y_prec;

converge:
VpChangeSign(y, 1);
#ifdef BIGDECIMAL_DEBUG
if (gfDebug) {
VpMult(r, y, y);
VpAddSub(f, x, r, -1);
printf("VpSqrt: iterations = %"PRIdSIZE"\n", nr);
VPrint(stdout, " y =% \n", y);
VPrint(stdout, " x =% \n", x);
VPrint(stdout, " x-y*y = % \n", f);
}
#endif /* BIGDECIMAL_DEBUG */
y->MaxPrec = y_prec;

Exit:
rbd_free_struct(f);
rbd_free_struct(r);
return 1;
}

/*
* Round relatively from the decimal point.
* f: rounding mode
Expand Down
3 changes: 0 additions & 3 deletions ext/bigdecimal/bigdecimal.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -195,7 +195,6 @@ typedef struct {
*/

#define VpBaseFig() BIGDECIMAL_COMPONENT_FIGURES
#define VpDblFig() BIGDECIMAL_DOUBLE_FIGURES

/* Zero,Inf,NaN (isinf(),isnan() used to check) */
VP_EXPORT double VpGetDoubleNaN(void);
Expand Down Expand Up @@ -229,8 +228,6 @@ VP_EXPORT void VpToString(Real *a, char *buf, size_t bufsize, size_t fFmt, int f
VP_EXPORT void VpToFString(Real *a, char *buf, size_t bufsize, size_t fFmt, int fPlus);
VP_EXPORT int VpCtoV(Real *a, const char *int_chr, size_t ni, const char *frac, size_t nf, const char *exp_chr, size_t ne);
VP_EXPORT int VpVtoD(double *d, SIGNED_VALUE *e, Real *m);
VP_EXPORT void VpDtoV(Real *m,double d);
VP_EXPORT int VpSqrt(Real *y,Real *x);
VP_EXPORT int VpActiveRound(Real *y, Real *x, unsigned short f, ssize_t il);
VP_EXPORT int VpMidRound(Real *y, unsigned short f, ssize_t nf);
VP_EXPORT int VpLeftRound(Real *y, unsigned short f, ssize_t nf);
Expand Down
33 changes: 33 additions & 0 deletions lib/bigdecimal.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,3 +3,36 @@
else
require 'bigdecimal.so'
end

class BigDecimal

# Returns the square root of the value.
#
# Result has at least prec significant digits.
#
def sqrt(prec)
if infinite? == 1
exception_mode = BigDecimal.mode(BigDecimal::EXCEPTION_ALL)
raise FloatDomainError, "Computation results in 'Infinity'" if exception_mode.anybits?(BigDecimal::EXCEPTION_INFINITY)
return INFINITY
end
raise ArgumentError, 'negative precision' if prec < 0
raise FloatDomainError, 'sqrt of negative value' if self < 0
raise FloatDomainError, "sqrt of 'NaN'(Not a Number)" if nan?

n_digits = n_significant_digits
prec = [prec, n_digits].max

@tompng tompng Jun 15, 2025

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This part is re-implementing #354


if n_digits < prec / 2
# Fast path for sqrt(16e100) => 4e50
ex = (n_digits - exponent + 1) / 2
n = (self * BigDecimal("1e#{2 * ex}")).to_i
sqrt = Integer.sqrt(n)
return BigDecimal(sqrt) * BigDecimal("1e#{-ex}") if sqrt * sqrt == n
end

ex = prec + BigDecimal.double_fig - exponent / 2
sqrt = Integer.sqrt(self * BigDecimal("1e#{2 * ex}"))
BigDecimal(sqrt) * BigDecimal("1e#{-ex}")
end
end
4 changes: 4 additions & 0 deletions test/bigdecimal/helper.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,10 @@ def teardown
end
end

def integer_sqrt_implemented?
RUBY_ENGINE != 'truffleruby' || Integer.sqrt(2**64 - 1) == 2**32 - 1
end

def under_gc_stress
stress, GC.stress = GC.stress, true
yield
Expand Down
28 changes: 26 additions & 2 deletions test/bigdecimal/test_bigdecimal.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1335,6 +1335,8 @@ def test_abs_bigdecimal
end

def test_sqrt_bigdecimal
pend 'Integer.sqrt not correctly implemented' unless integer_sqrt_implemented?

x = BigDecimal("0.09")
assert_in_delta(0.3, x.sqrt(1), 0.001)
x = BigDecimal((2**100).to_s)
Expand All @@ -1345,8 +1347,6 @@ def test_sqrt_bigdecimal
assert_equal(true, (x.sqrt(300) - y).abs < BigDecimal("1E#{e-300}"))
x = BigDecimal("-" + (2**100).to_s)
assert_raise_with_message(FloatDomainError, "sqrt of negative value") { x.sqrt(1) }
x = BigDecimal((2**200).to_s)
assert_equal(2**100, x.sqrt(1))

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

Expand All @@ -1359,9 +1359,17 @@ def test_sqrt_bigdecimal
assert_equal(0, BigDecimal("-0").sqrt(1))
assert_equal(1, BigDecimal("1").sqrt(1))
assert_positive_infinite(BigDecimal("Infinity").sqrt(1))

assert_equal(BigDecimal('11.1'), BigDecimal('123.21').sqrt(100))
assert_equal(BigDecimal('11e20'), BigDecimal('121e40').sqrt(100))
assert_in_epsilon(Math.sqrt(121e41), BigDecimal('121e41').sqrt(100))
assert_in_epsilon(Math.sqrt(121.5e40), BigDecimal('121.5e40').sqrt(100))
assert_in_epsilon(Math.sqrt(2e100), BigDecimal('2e100').sqrt(10))
end

def test_sqrt_5266
pend 'Integer.sqrt not correctly implemented' unless integer_sqrt_implemented?

x = BigDecimal('2' + '0'*100)
assert_equal('0.14142135623730950488016887242096980785696718753769480731',
x.sqrt(56).to_s(56).split(' ')[0])
Expand All @@ -1375,6 +1383,22 @@ def test_sqrt_5266
x.sqrt(109).to_s(109).split(' ')[0])
end

def test_sqrt_minimum_precision
pend 'Integer.sqrt not correctly implemented' unless integer_sqrt_implemented?

x = BigDecimal((2**200).to_s)
assert_equal(2**100, x.sqrt(1))

x = BigDecimal('1' * 60 + '.' + '1' * 40)
assert_in_delta(BigDecimal('3' * 30 + '.' + '3' * 70), x.sqrt(1), BigDecimal('1e-70'))

x = BigDecimal('1' * 40 + '.' + '1' * 60)
assert_in_delta(BigDecimal('3' * 20 + '.' + '3' * 80), x.sqrt(1), BigDecimal('1e-80'))

x = BigDecimal('0.' + '0' * 50 + '1' * 100)
assert_in_delta(BigDecimal('0.' + '0' * 25 + '3' * 100), x.sqrt(1), BigDecimal('1e-125'))
end

def test_fix
x = BigDecimal("1.1")
assert_equal(1, x.fix)
Expand Down
4 changes: 4 additions & 0 deletions test/bigdecimal/test_bigmath.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@ def test_const
end

def test_sqrt
pend 'Integer.sqrt not correctly implemented' unless integer_sqrt_implemented?

@tompng tompng Jun 18, 2025

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pend test_sqrt and test_atan that uses Integer.sqrt internally. Workaround for ci failure in truffleruby.
Integer.sqrt is fixed in truffleruby-head.

Another options are:

  • Override Integer.sqrt in test
  • Implement _integer_sqrt(n) = Integer.sqrt(n) or ruby_implementation in bigdecimal


assert_in_delta(2**0.5, sqrt(BigDecimal("2"), N))
assert_equal(10, sqrt(BigDecimal("100"), N))
assert_equal(0.0, sqrt(BigDecimal("0"), N))
Expand Down Expand Up @@ -81,6 +83,8 @@ def test_cos
end

def test_atan
pend 'Integer.sqrt not correctly implemented' unless integer_sqrt_implemented?

assert_equal(0.0, atan(BigDecimal("0.0"), N))
assert_in_delta(Math::PI/4, atan(BigDecimal("1.0"), N))
assert_in_delta(Math::PI/6, atan(sqrt(BigDecimal("3.0"), N) / 3, N))
Expand Down
Loading