diff options
| -rw-r--r-- | bignum.c | 216 | ||||
| -rw-r--r-- | internal.h | 54 | ||||
| -rw-r--r-- | numeric.c | 78 | ||||
| -rw-r--r-- | test/ruby/test_numeric.rb | 14 |
4 files changed, 302 insertions, 60 deletions
@@ -6877,6 +6877,222 @@ rb_big_isqrt(VALUE n) return x; } +#ifdef USE_GMP +static void +bary_powm_gmp(BDIGIT *zds, size_t zn, const BDIGIT *xds, size_t xn, const BDIGIT *yds, size_t yn, const BDIGIT *mds, size_t mn) +{ + const size_t nails = (sizeof(BDIGIT)-SIZEOF_BDIGIT)*CHAR_BIT; + mpz_t z, x, y, m; + size_t count; + mpz_init(x); + mpz_init(y); + mpz_init(m); + mpz_init(z); + mpz_import(x, xn, -1, sizeof(BDIGIT), 0, nails, xds); + mpz_import(y, yn, -1, sizeof(BDIGIT), 0, nails, yds); + mpz_import(m, mn, -1, sizeof(BDIGIT), 0, nails, mds); + mpz_powm(z, x, y, m); + mpz_export(zds, &count, -1, sizeof(BDIGIT), 0, nails, z); + BDIGITS_ZERO(zds+count, zn-count); + mpz_clear(x); + mpz_clear(y); + mpz_clear(m); + mpz_clear(z); +} +#endif + +static VALUE +int_pow_tmp3(VALUE x, VALUE y, VALUE m, int nega_flg) +{ +#ifdef USE_GMP + VALUE z; + size_t xn, yn, mn, zn; + + if (FIXNUM_P(x)) { + x = rb_int2big(FIX2LONG(x)); + } + if (FIXNUM_P(y)) { + y = rb_int2big(FIX2LONG(y)); + } + assert(RB_BIGNUM_TYPE_P(m)); + xn = BIGNUM_LEN(x); + yn = BIGNUM_LEN(y); + mn = BIGNUM_LEN(m); + zn = mn; + z = bignew(zn, 1); + bary_powm_gmp(BDIGITS(z), zn, BDIGITS(x), xn, BDIGITS(y), yn, BDIGITS(m), mn); + if (nega_flg & BIGNUM_POSITIVE_P(z)) { + z = rb_funcall(z, '-', 1, m); + } + RB_GC_GUARD(x); + RB_GC_GUARD(y); + RB_GC_GUARD(m); + return rb_big_norm(z); +#else + VALUE tmp = LONG2FIX(1L); + long yy; + + for (/*NOP*/; ! FIXNUM_P(y); y = rb_funcall(y, rb_intern(">>"), 1, LONG2FIX(1L))) { + if (RTEST(rb_funcall(y, rb_intern("odd?"), 0))) { + tmp = rb_funcall(tmp, '*', 1, x); + tmp = rb_int_modulo(tmp, m); + } + x = rb_funcall(x, '*', 1, x); + x = rb_int_modulo(x, m); + } + for (yy = FIX2LONG(y); yy; yy >>= 1L) { + if (yy & 1L) { + tmp = rb_funcall(tmp, '*', 1, x); + tmp = rb_int_modulo(tmp, m); + } + x = rb_funcall(x, '*', 1, x); + x = rb_int_modulo(x, m); + } + + if (nega_flg && RTEST(rb_funcall(tmp, rb_intern("positive?"), 0))) { + tmp = rb_funcall(tmp, '-', 1, m); + } + return tmp; +#endif +} + +/* + * Integer#pow + */ + +static VALUE +int_pow_tmp1(VALUE x, VALUE y, long mm, int nega_flg) +{ + long xx = FIX2LONG(x); + long tmp = 1L; + long yy; + + for (/*NOP*/; ! FIXNUM_P(y); y = rb_funcall(y, idGTGT, 1, LONG2FIX(1L))) { + if (RTEST(rb_int_odd_p(y))) { + tmp = (tmp * xx) % mm; + } + xx = (xx * xx) % mm; + } + for (yy = FIX2LONG(y); yy; yy >>= 1L) { + if (yy & 1L) { + tmp = (tmp * xx) % mm; + } + xx = (xx * xx) % mm; + } + + if (nega_flg && tmp) { + tmp -= mm; + } + return LONG2FIX(tmp); +} + +static VALUE +int_pow_tmp2(VALUE x, VALUE y, long mm, int nega_flg) +{ + long tmp = 1L; + long yy; +#ifdef DLONG + DLONG const mmm = mm; + long xx = FIX2LONG(x); + + for (/*NOP*/; ! FIXNUM_P(y); y = rb_funcall(y, idGTGT, 1, LONG2FIX(1L))) { + if (RTEST(rb_int_odd_p(y))) { + tmp = ((DLONG)tmp * (DLONG)xx) % mmm; + } + xx = ((DLONG)xx * (DLONG)xx) % mmm; + } + for (yy = FIX2LONG(y); yy; yy >>= 1L) { + if (yy & 1L) { + tmp = ((DLONG)tmp * (DLONG)xx) % mmm; + } + xx = ((DLONG)xx * (DLONG)xx) % mmm; + } +#else + VALUE const m = LONG2FIX(mm); + VALUE tmp2 = LONG2FIX(tmp); + + for (/*NOP*/; ! FIXNUM_P(y); y = rb_funcall(y, idGTGT, 1, LONG2FIX(1L))) { + if (RTEST(rb_int_odd_p(y))) { + tmp2 = rb_fix_mul_fix(tmp2, x); + tmp2 = rb_int_modulo(tmp2, m); + } + x = rb_fix_mul_fix(x, x); + x = rb_int_modulo(x, m); + } + for (yy = FIX2LONG(y); yy; yy >>= 1L) { + if (yy & 1L) { + tmp2 = rb_fix_mul_fix(tmp2, x); + tmp2 = rb_int_modulo(tmp2, m); + } + x = rb_fix_mul_fix(x, x); + x = rb_int_modulo(x, m); + } + + tmp = FIX2LONG(tmp2); +#endif + if (nega_flg && tmp) { + tmp -= mm; + } + return LONG2FIX(tmp); +} + +/* + * Document-method: Integer#pow + * call-seq: + * integer.pow(integer) -> integer + * integer.pow(integer, integer) -> integer + * + * Returns (modular) exponentiation as: + * + * a.pow(b) #=> same as a**b + * a.pow(b, m) #=> same as (a**b) % m, but doesn't make huge values as temporary + */ +VALUE +rb_int_powm(int const argc, VALUE * const argv, VALUE const num) +{ + rb_check_arity(argc, 1, 2); + + if (argc == 1) { + return rb_funcall(num, rb_intern("**"), 1, argv[0]); + } + else { + VALUE const a = num; + VALUE const b = argv[0]; + VALUE m = argv[1]; + int nega_flg = 0; + if ( ! RB_INTEGER_TYPE_P(b)) { + rb_raise(rb_eTypeError, "Integer#pow() 2nd argument not allowed unless a 1st argument is integer"); + } + if (rb_num_negative_int_p(b)) { + rb_raise(rb_eRangeError, "Integer#pow() 1st argument cannot be negative when 2nd argument specified"); + } + if (!RB_INTEGER_TYPE_P(m)) { + rb_raise(rb_eTypeError, "Integer#pow() 2nd argument not allowed unless all arguments are integers"); + } + + if (rb_num_negative_int_p(m)) { + m = rb_funcall(m, idUMinus, 0); + nega_flg = 1; + } + + if (!rb_num_positive_int_p(m)) { + rb_num_zerodiv(); + } + if (FIXNUM_P(m)) { + long const half_val = (long)HALF_LONG_MSB; + long const mm = FIX2LONG(m); + if (mm <= half_val) { + return int_pow_tmp1(rb_int_modulo(a, m), b, mm, nega_flg); + } else { + return int_pow_tmp2(rb_int_modulo(a, m), b, mm, nega_flg); + } + } else if (RB_TYPE_P(m, T_BIGNUM)) { + return int_pow_tmp3(rb_int_modulo(a, m), b, m, nega_flg); + } + } + UNREACHABLE; +} + /* * Bignum objects hold integers outside the range of * Fixnum. Bignum objects are created diff --git a/internal.h b/internal.h index 9d33ed99d1..81ecb876ad 100644 --- a/internal.h +++ b/internal.h @@ -39,6 +39,12 @@ extern "C" { # endif #endif +/* The most significant bit of the lower part of half-long integer. + * If sizeof(long) == 4, this is 0x8000. + * If sizeof(long) == 8, this is 0x80000000. + */ +#define HALF_LONG_MSB ((SIGNED_VALUE)1<<((SIZEOF_LONG*CHAR_BIT-1)/2)) + #define LIKELY(x) RB_LIKELY(x) #define UNLIKELY(x) RB_UNLIKELY(x) @@ -1074,6 +1080,7 @@ VALUE rb_big_gt(VALUE x, VALUE y); VALUE rb_big_ge(VALUE x, VALUE y); VALUE rb_big_lt(VALUE x, VALUE y); VALUE rb_big_le(VALUE x, VALUE y); +VALUE rb_int_powm(int const argc, VALUE * const argv, VALUE const num); /* class.c */ VALUE rb_class_boot(VALUE); @@ -1380,6 +1387,53 @@ VALUE rb_int_and(VALUE x, VALUE y); VALUE rb_int_lshift(VALUE x, VALUE y); VALUE rb_int_div(VALUE x, VALUE y); VALUE rb_int_abs(VALUE num); +VALUE rb_int_odd_p(VALUE num); + +static inline VALUE +rb_num_compare_with_zero(VALUE num, ID mid) +{ + VALUE zero = INT2FIX(0); + VALUE r = rb_check_funcall(num, mid, 1, &zero); + if (r == Qundef) { + rb_cmperr(num, zero); + } + return r; +} + +static inline int +rb_num_positive_int_p(VALUE num) +{ + const ID mid = '>'; + + if (FIXNUM_P(num)) { + if (rb_method_basic_definition_p(rb_cInteger, mid)) + return FIXNUM_POSITIVE_P(num); + } + else if (RB_TYPE_P(num, T_BIGNUM)) { + if (rb_method_basic_definition_p(rb_cInteger, mid)) + return BIGNUM_POSITIVE_P(num); + } + return RTEST(rb_num_compare_with_zero(num, mid)); +} + + +static inline int +rb_num_negative_int_p(VALUE num) +{ + const ID mid = '<'; + + if (FIXNUM_P(num)) { + if (rb_method_basic_definition_p(rb_cInteger, mid)) + return FIXNUM_NEGATIVE_P(num); + } + else if (RB_TYPE_P(num, T_BIGNUM)) { + if (rb_method_basic_definition_p(rb_cInteger, mid)) + return BIGNUM_NEGATIVE_P(num); + } + return RTEST(rb_num_compare_with_zero(num, mid)); +} + + VALUE rb_float_abs(VALUE flt); VALUE rb_float_equal(VALUE x, VALUE y); VALUE rb_float_eql(VALUE x, VALUE y); @@ -162,7 +162,6 @@ static VALUE fix_mul(VALUE x, VALUE y); static VALUE fix_lshift(long, unsigned long); static VALUE fix_rshift(long, unsigned long); static VALUE int_pow(long x, unsigned long y); -static VALUE int_odd_p(VALUE x); static VALUE int_even_p(VALUE x); static int int_round_zero_p(VALUE num, int ndigits); VALUE rb_int_floor(VALUE num, int ndigits); @@ -271,17 +270,6 @@ rb_num_to_uint(VALUE val, unsigned int *ret) #define method_basic_p(klass) rb_method_basic_definition_p(klass, mid) -static VALUE -compare_with_zero(VALUE num, ID mid) -{ - VALUE zero = INT2FIX(0); - VALUE r = rb_check_funcall(num, mid, 1, &zero); - if (r == Qundef) { - rb_cmperr(num, zero); - } - return r; -} - static inline int int_pos_p(VALUE num) { @@ -306,42 +294,10 @@ int_neg_p(VALUE num) rb_raise(rb_eTypeError, "not an Integer"); } -static inline int -positive_int_p(VALUE num) -{ - const ID mid = '>'; - - if (FIXNUM_P(num)) { - if (method_basic_p(rb_cInteger)) - return FIXNUM_POSITIVE_P(num); - } - else if (RB_TYPE_P(num, T_BIGNUM)) { - if (method_basic_p(rb_cInteger)) - return BIGNUM_POSITIVE_P(num); - } - return RTEST(compare_with_zero(num, mid)); -} - -static inline int -negative_int_p(VALUE num) -{ - const ID mid = '<'; - - if (FIXNUM_P(num)) { - if (method_basic_p(rb_cInteger)) - return FIXNUM_NEGATIVE_P(num); - } - else if (RB_TYPE_P(num, T_BIGNUM)) { - if (method_basic_p(rb_cInteger)) - return BIGNUM_NEGATIVE_P(num); - } - return RTEST(compare_with_zero(num, mid)); -} - int rb_num_negative_p(VALUE num) { - return negative_int_p(num); + return rb_num_negative_int_p(num); } static VALUE @@ -665,10 +621,10 @@ num_remainder(VALUE x, VALUE y) VALUE z = num_funcall1(x, '%', y); if ((!rb_equal(z, INT2FIX(0))) && - ((negative_int_p(x) && - positive_int_p(y)) || - (positive_int_p(x) && - negative_int_p(y)))) { + ((rb_num_negative_int_p(x) && + rb_num_positive_int_p(y)) || + (rb_num_positive_int_p(x) && + rb_num_negative_int_p(y)))) { return rb_funcall(z, '-', 1, y); } return z; @@ -769,7 +725,7 @@ num_int_p(VALUE num) static VALUE num_abs(VALUE num) { - if (negative_int_p(num)) { + if (rb_num_negative_int_p(num)) { return num_funcall0(num, idUMinus); } return num; @@ -886,7 +842,7 @@ num_positive_p(VALUE num) if (method_basic_p(rb_cInteger)) return BIGNUM_POSITIVE_P(num) && !rb_bigzero_p(num) ? Qtrue : Qfalse; } - return compare_with_zero(num, mid); + return rb_num_compare_with_zero(num, mid); } /* @@ -899,7 +855,7 @@ num_positive_p(VALUE num) static VALUE num_negative_p(VALUE num) { - return negative_int_p(num) ? Qtrue : Qfalse; + return rb_num_negative_int_p(num) ? Qtrue : Qfalse; } @@ -2072,7 +2028,7 @@ int_round_half_down(SIGNED_VALUE x, SIGNED_VALUE y) static int int_half_p_half_even(VALUE num, VALUE n, VALUE f) { - return (int)int_odd_p(rb_int_idiv(n, f)); + return (int)rb_int_odd_p(rb_int_idiv(n, f)); } static int @@ -2939,7 +2895,7 @@ rb_fix2uint(VALUE val) } num = FIX2ULONG(val); - check_uint(num, negative_int_p(val)); + check_uint(num, rb_num_negative_int_p(val)); return num; } #else @@ -3025,7 +2981,7 @@ rb_fix2ushort(VALUE val) } num = FIX2ULONG(val); - check_ushort(num, negative_int_p(val)); + check_ushort(num, rb_num_negative_int_p(val)); return num; } @@ -3168,8 +3124,8 @@ int_int_p(VALUE num) * Returns +true+ if +int+ is an odd number. */ -static VALUE -int_odd_p(VALUE num) +VALUE +rb_int_odd_p(VALUE num) { if (FIXNUM_P(num)) { if (num & 2) { @@ -3567,7 +3523,7 @@ rb_int_minus(VALUE x, VALUE y) } -#define SQRT_LONG_MAX ((SIGNED_VALUE)1<<((SIZEOF_LONG*CHAR_BIT-1)/2)) +#define SQRT_LONG_MAX HALF_LONG_MSB /*tests if N*N would overflow*/ #define FIT_SQRT_LONG(n) (((n)<SQRT_LONG_MAX)&&((n)>=-SQRT_LONG_MAX)) @@ -3976,7 +3932,7 @@ fix_pow(VALUE x, VALUE y) if (int_even_p(y)) return INT2FIX(1); else return INT2FIX(-1); } - if (negative_int_p(y)) + if (rb_num_negative_int_p(y)) return num_funcall1(rb_rational_raw1(x), idPow, y); if (a == 0) return INT2FIX(0); x = rb_int2big(FIX2LONG(x)); @@ -5394,7 +5350,7 @@ Init_Numeric(void) rb_define_method(rb_cInteger, "to_s", int_to_s, -1); rb_define_alias(rb_cInteger, "inspect", "to_s"); rb_define_method(rb_cInteger, "integer?", int_int_p, 0); - rb_define_method(rb_cInteger, "odd?", int_odd_p, 0); + rb_define_method(rb_cInteger, "odd?", rb_int_odd_p, 0); rb_define_method(rb_cInteger, "even?", int_even_p, 0); rb_define_method(rb_cInteger, "upto", int_upto, 1); rb_define_method(rb_cInteger, "downto", int_downto, 1); @@ -5426,6 +5382,8 @@ Init_Numeric(void) rb_define_method(rb_cInteger, "fdiv", rb_int_fdiv, 1); rb_define_method(rb_cInteger, "**", rb_int_pow, 1); + rb_define_method(rb_cInteger, "pow", rb_int_powm, -1); + rb_define_method(rb_cInteger, "abs", rb_int_abs, 0); rb_define_method(rb_cInteger, "magnitude", rb_int_abs, 0); diff --git a/test/ruby/test_numeric.rb b/test/ruby/test_numeric.rb index 49008247ed..9b17fbb2a5 100644 --- a/test/ruby/test_numeric.rb +++ b/test/ruby/test_numeric.rb @@ -384,4 +384,18 @@ class TestNumeric < Test::Unit::TestCase end end end + + def test_pow + assert_equal(3**3 % 8, 3.pow(3, 8)) + assert_equal(3**3 % -8, 3.pow(3,-8)) + assert_equal(3**2 % -2, 3.pow(2,-2)) + assert_equal((-3)**3 % 8, -3.pow(3,8)) + assert_equal((-3)**3 % -8, -3.pow(3,-8)) + assert_equal(5**2 % -8, 5.pow(2,-8)) + assert_equal(4481650795473624846969600733813414725093, + 2120078484650058507891187874713297895455. + pow(5478118174010360425845660566650432540723, + 5263488859030795548286226023720904036518)) + end + end |