Skip to content

NTT multiplication and Newton-Raphson division #407

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

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

NTT multiplication and Newton-Raphson division #407

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
61820c4
Implement faster multiplication using Number Theoretic Transform
tompng Aug 17, 2025
2513f01
Implement Newton-Raphson division
tompng Aug 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
2 changes: 2 additions & 0 deletions bigdecimal.gemspec
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,9 +43,11 @@ Gem::Specification.new do |s|
ext/bigdecimal/bigdecimal.c
ext/bigdecimal/bigdecimal.h
ext/bigdecimal/bits.h
ext/bigdecimal/div.h
ext/bigdecimal/feature.h
ext/bigdecimal/missing.c
ext/bigdecimal/missing.h
ext/bigdecimal/ntt.h
ext/bigdecimal/missing/dtoa.c
ext/bigdecimal/static_assert.h
]
Expand Down
86 changes: 74 additions & 12 deletions ext/bigdecimal/bigdecimal.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,10 +29,18 @@
#endif

#include "bits.h"
#include "div.h"
#include "static_assert.h"

#define BIGDECIMAL_VERSION "3.3.0"

#if SIZEOF_DECDIG == 4
#define USE_NTT_MULTIPLICATION 1
#include "ntt.h"
#define NTT_MULTIPLICATION_THRESHOLD 100
#define NEWTON_RAPHSON_DIVISION_THRESHOLD 200
#endif

/* #define ENABLE_NUMERIC_STRING */

#define SIGNED_VALUE_MAX INTPTR_MAX
Expand Down Expand Up @@ -75,11 +83,6 @@ static struct {
uint8_t mode;
} rbd_rounding_modes[RBD_NUM_ROUNDING_MODES];

typedef struct {
VALUE bigdecimal;
Real *real;
} BDVALUE;

typedef struct {
VALUE bigdecimal_or_nil;
Real *real_or_null;
Expand Down Expand Up @@ -207,7 +210,6 @@ rbd_allocate_struct_zero(int sign, size_t const digits)
static unsigned short VpGetException(void);
static void VpSetException(unsigned short f);
static void VpCheckException(Real *p, bool always);
static int AddExponent(Real *a, SIGNED_VALUE n);
static VALUE CheckGetValue(BDVALUE v);
static void VpInternalRound(Real *c, size_t ixDigit, DECDIG vPrev, DECDIG v);
static int VpLimitRound(Real *c, size_t ixDigit);
Expand Down Expand Up @@ -1112,9 +1114,6 @@ BigDecimal_check_num(Real *p)
VpCheckException(p, true);
}

static VALUE BigDecimal_fix(VALUE self);
static VALUE BigDecimal_split(VALUE self);

/* Returns the value as an Integer.
*
* If the BigDecimal is infinity or NaN, raises FloatDomainError.
Expand Down Expand Up @@ -3256,19 +3255,39 @@ BigDecimal_literal(const char *str)

#ifdef BIGDECIMAL_USE_VP_TEST_METHODS
VALUE
BigDecimal_vpdivd(VALUE self, VALUE r, VALUE cprec) {
BDVALUE a,b,c,d;
BigDecimal_vpdivd_generic(VALUE self, VALUE r, VALUE cprec, void (*vpdivd_func)(Real*, Real*, Real*, Real*)) {
BDVALUE a, b, c, d;
size_t cn = NUM2INT(cprec);
a = GetBDValueMust(self);
b = GetBDValueMust(r);
c = NewZeroWrap(1, cn * BASE_FIG);
d = NewZeroWrap(1, VPDIVD_REM_PREC(a.real, b.real, c.real) * BASE_FIG);
VpDivd(c.real, d.real, a.real, b.real);
vpdivd_func(c.real, d.real, a.real, b.real);
RB_GC_GUARD(a.bigdecimal);
RB_GC_GUARD(b.bigdecimal);
return rb_assoc_new(c.bigdecimal, d.bigdecimal);
}

void
VpDivdNormal(Real *c, Real *r, Real *a, Real *b) {
VpDivd(c, r, a, b);
}

VALUE
BigDecimal_vpdivd(VALUE self, VALUE r, VALUE cprec) {
return BigDecimal_vpdivd_generic(self, r, cprec, VpDivdNormal);
}

VALUE
BigDecimal_vpdivd_newton(VALUE self, VALUE r, VALUE cprec) {
return BigDecimal_vpdivd_generic(self, r, cprec, VpDivdNewton);
}

VALUE
BigDecimal_newton_raphson_inverse(VALUE self, VALUE prec) {
return newton_raphson_inverse(self, NUM2SIZET(prec));
}

VALUE
BigDecimal_vpmult(VALUE self, VALUE v) {
BDVALUE a,b,c;
Expand All @@ -3280,6 +3299,25 @@ BigDecimal_vpmult(VALUE self, VALUE v) {
RB_GC_GUARD(b.bigdecimal);
return c.bigdecimal;
}

#if SIZEOF_DECDIG == 4
VALUE
BigDecimal_nttmult(VALUE self, VALUE v) {
BDVALUE a,b,c;
a = GetBDValueMust(self);
b = GetBDValueMust(v);
c = NewZeroWrap(1, VPMULT_RESULT_PREC(a.real, b.real) * BASE_FIG);
ntt_multiply(a.real->Prec, b.real->Prec, a.real->frac, b.real->frac, c.real->frac);
VpSetSign(c.real, a.real->sign * b.real->sign);
c.real->exponent = a.real->exponent + b.real->exponent;
c.real->Prec = a.real->Prec + b.real->Prec;
VpNmlz(c.real);
RB_GC_GUARD(a.bigdecimal);
RB_GC_GUARD(b.bigdecimal);
return c.bigdecimal;
}
#endif

#endif /* BIGDECIMAL_USE_VP_TEST_METHODS */

/* Document-class: BigDecimal
Expand Down Expand Up @@ -3651,7 +3689,12 @@ Init_bigdecimal(void)

#ifdef BIGDECIMAL_USE_VP_TEST_METHODS
rb_define_method(rb_cBigDecimal, "vpdivd", BigDecimal_vpdivd, 2);
rb_define_method(rb_cBigDecimal, "vpdivd_newton", BigDecimal_vpdivd_newton, 2);
rb_define_method(rb_cBigDecimal, "newton_raphson_inverse", BigDecimal_newton_raphson_inverse, 1);
rb_define_method(rb_cBigDecimal, "vpmult", BigDecimal_vpmult, 1);
#ifdef USE_NTT_MULTIPLICATION
rb_define_method(rb_cBigDecimal, "nttmult", BigDecimal_nttmult, 1);
#endif
#endif /* BIGDECIMAL_USE_VP_TEST_METHODS */

#define ROUNDING_MODE(i, name, value) \
Expand Down Expand Up @@ -4934,6 +4977,15 @@ VpMult(Real *c, Real *a, Real *b)
c->exponent = a->exponent; /* set exponent */
VpSetSign(c, VpGetSign(a) * VpGetSign(b)); /* set sign */
if (!AddExponent(c, b->exponent)) return 0;

#ifdef USE_NTT_MULTIPLICATION
if (b->Prec >= NTT_MULTIPLICATION_THRESHOLD) {
ntt_multiply((uint32_t)a->Prec, (uint32_t)b->Prec, a->frac, b->frac, c->frac);
c->Prec = a->Prec + b->Prec;
goto Cleanup;
}
#endif

carry = 0;
nc = ind_c = MxIndAB;
memset(c->frac, 0, (nc + 1) * sizeof(DECDIG)); /* Initialize c */
Expand Down Expand Up @@ -4980,6 +5032,8 @@ VpMult(Real *c, Real *a, Real *b)
}
}
}

Cleanup:
VpNmlz(c);

Exit:
Expand Down Expand Up @@ -5027,6 +5081,14 @@ VpDivd(Real *c, Real *r, Real *a, Real *b)

if (word_a > word_r || word_b + word_c - 2 >= word_r) goto space_error;

#ifdef USE_NTT_MULTIPLICATION
// Newton-Raphson division requires multiplication to be faster than O(n^2)
if (word_c >= NEWTON_RAPHSON_DIVISION_THRESHOLD && word_b >= NEWTON_RAPHSON_DIVISION_THRESHOLD) {
VpDivdNewton(c, r, a, b);
goto Exit;
}
#endif

for (i = 0; i < word_a; ++i) r->frac[i] = a->frac[i];
for (i = word_a; i < word_r; ++i) r->frac[i] = 0;
for (i = 0; i < word_c; ++i) c->frac[i] = 0;
Expand Down
26 changes: 26 additions & 0 deletions ext/bigdecimal/bigdecimal.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -188,6 +188,11 @@ typedef struct {
DECDIG frac[FLEXIBLE_ARRAY_SIZE]; /* Array of fraction part. */
} Real;

typedef struct {
VALUE bigdecimal;
Real *real;
} BDVALUE;

/*
* ------------------
* EXPORTables.
Expand Down Expand Up @@ -232,10 +237,31 @@ 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);
VP_EXPORT void VpFrac(Real *y, Real *x);
VP_EXPORT int AddExponent(Real *a, SIGNED_VALUE n);

/* VP constants */
VP_EXPORT Real *VpOne(void);

/*
* **** BigDecimal part ****
*/
VP_EXPORT VALUE BigDecimal_lt(VALUE self, VALUE r);
VP_EXPORT VALUE BigDecimal_ge(VALUE self, VALUE r);
VP_EXPORT VALUE BigDecimal_exponent(VALUE self);
VP_EXPORT VALUE BigDecimal_fix(VALUE self);
VP_EXPORT VALUE BigDecimal_frac(VALUE self);
VP_EXPORT VALUE BigDecimal_add(VALUE self, VALUE b);
VP_EXPORT VALUE BigDecimal_sub(VALUE self, VALUE b);
VP_EXPORT VALUE BigDecimal_mult(VALUE self, VALUE b);
VP_EXPORT VALUE BigDecimal_add2(VALUE self, VALUE b, VALUE n);
VP_EXPORT VALUE BigDecimal_sub2(VALUE self, VALUE b, VALUE n);
VP_EXPORT VALUE BigDecimal_mult2(VALUE self, VALUE b, VALUE n);
VP_EXPORT VALUE BigDecimal_split(VALUE self);
VP_EXPORT VALUE BigDecimal_decimal_shift(VALUE self, VALUE v);
VP_EXPORT inline BDVALUE GetBDValueMust(VALUE v);
VP_EXPORT inline BDVALUE rbd_allocate_struct_zero_wrap(int sign, size_t const digits);
#define NewZeroWrap rbd_allocate_struct_zero_wrap

/*
* ------------------
* MACRO definitions.
Expand Down
Loading