1berry/ruby
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

[Bug #11183] Fix rb_complex_pow for special angles

Browse Source 
Add a special treatment for when the argument of self is an
integral multiple of 45 degrees.

  1i ** (10 ** 100)         #=> 1+0i
  1i ** (10 ** 100 + 1)     #=> 0+1i
  (1+1i) ** (10 ** 100)     # warning: in a**b, b may be too big
                            #=> (Infinity+0.0i)
  (1+1i) ** (10 ** 100 + 1) # warning: in a**b, b may be too big
                            #=> (Infinity+Infinity*i)
This commit is contained in:
Kouhei Yanagita 2023-11-16 20:32:53 +09:00 committed by Nobuyoshi Nakada
parent b6b31f673d
commit 04eb40b633
2 changed files with 151 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

86
complex.c
View File

@ -982,6 +982,84 @@ f_reciprocal(VALUE x)
return f_quo(ONE, x);
}
static VALUE
zero_for(VALUE x)
{
if (RB_FLOAT_TYPE_P(x))
return DBL2NUM(0);
if (RB_TYPE_P(x, T_RATIONAL))
return rb_rational_new(INT2FIX(0), INT2FIX(1));
return INT2FIX(0);
}
static VALUE
complex_pow_for_special_angle(VALUE self, VALUE other)
{
if (!rb_integer_type_p(other)) {
return Qundef;
}
get_dat1(self);
VALUE x = Qundef;
int dir;
if (f_zero_p(dat->imag)) {
x = dat->real;
dir = 0;
}
else if (f_zero_p(dat->real)) {
x = dat->imag;
dir = 2;
}
else if (f_eqeq_p(dat->real, dat->imag)) {
x = dat->real;
dir = 1;
}
else if (f_eqeq_p(dat->real, f_negate(dat->imag))) {
x = dat->imag;
dir = 3;
}
if (x == Qundef) return x;
if (f_negative_p(x)) {
x = f_negate(x);
dir += 4;
}
VALUE zx;
if (dir % 2 == 0) {
zx = rb_num_pow(x, other);
}
else {
zx = rb_num_pow(
rb_funcall(rb_int_mul(TWO, x), '*', 1, x),
rb_int_div(other, TWO)
);
if (rb_int_odd_p(other)) {
zx = rb_funcall(zx, '*', 1, x);
}
}
static const int dirs[][2] = {
{1, 0}, {1, 1}, {0, 1}, {-1, 1}, {-1, 0}, {-1, -1}, {0, -1}, {1, -1}
};
int z_dir = FIX2INT(rb_int_modulo(rb_int_mul(INT2FIX(dir), other), INT2FIX(8)));
VALUE zr, zi;
switch (dirs[z_dir][0]) {
case 0: zr = zero_for(zx); break;
case 1: zr = zx; break;
case -1: zr = f_negate(zx); break;
}
switch (dirs[z_dir][1]) {
case 0: zi = zero_for(zx); break;
case 1: zi = zx; break;
case -1: zi = f_negate(zx); break;
}
return nucomp_s_new_internal(CLASS_OF(self), zr, zi);
}
/*
* call-seq:
* cmp ** numeric -> complex
@ -1007,6 +1085,14 @@ rb_complex_pow(VALUE self, VALUE other)
other = dat->real; /* c14n */
}
if (other == ONE) {
get_dat1(self);
return nucomp_s_new_internal(CLASS_OF(self), dat->real, dat->imag);
}
VALUE result = complex_pow_for_special_angle(self, other);
if (result != Qundef) return result;
if (RB_TYPE_P(other, T_COMPLEX)) {
VALUE r, theta, nr, ntheta;

65
test/ruby/test_complex.rb
View File

@ -526,6 +526,71 @@ class Complex_Test < Test::Unit::TestCase
r = c ** Rational(-2,3)
assert_in_delta(0.432, r.real, 0.001)
assert_in_delta(-0.393, r.imag, 0.001)
end
def test_expt_for_special_angle
c = Complex(1, 0) ** 100000000000000000000000000000000
assert_equal(Complex(1, 0), c)
c = Complex(-1, 0) ** 10000000000000000000000000000000
assert_equal(Complex(1, 0), c)
c = Complex(-1, 0) ** 10000000000000000000000000000001
assert_equal(Complex(-1, 0), c)
c = Complex(0, 1) ** 100000000000000000000000000000000
assert_equal(Complex(1, 0), c)
c = Complex(0, 1) ** 100000000000000000000000000000001
assert_equal(Complex(0, 1), c)
c = Complex(0, 1) ** 100000000000000000000000000000002
assert_equal(Complex(-1, 0), c)
c = Complex(0, 1) ** 100000000000000000000000000000003
assert_equal(Complex(0, -1), c)
c = Complex(0, -1) ** 100000000000000000000000000000000
assert_equal(Complex(1, 0), c)
c = Complex(0, -1) ** 100000000000000000000000000000001
assert_equal(Complex(0, -1), c)
c = Complex(0, -1) ** 100000000000000000000000000000002
assert_equal(Complex(-1, 0), c)
c = Complex(0, -1) ** 100000000000000000000000000000003
assert_equal(Complex(0, 1), c)
c = Complex(1, 1) ** 1
assert_equal(Complex(1, 1), c)
c = Complex(1, 1) ** 2
assert_equal(Complex(0, 2), c)
c = Complex(1, 1) ** 3
assert_equal(Complex(-2, 2), c)
c = Complex(1, 1) ** 4
assert_equal(Complex(-4, 0), c)
c = Complex(1, 1) ** 5
assert_equal(Complex(-4, -4), c)
c = Complex(1, 1) ** 6
assert_equal(Complex(0, -8), c)
c = Complex(1, 1) ** 7
assert_equal(Complex(8, -8), c)
c = Complex(-2, -2) ** 3
assert_equal(Complex(16, -16), c)
c = Complex(2, -2) ** 3
assert_equal(Complex(-16, -16), c)
c = Complex(-2, 2) ** 3
assert_equal(Complex(16, 16), c)
c = Complex(0.0, -888888888888888.0)**8888
assert_not_predicate(c.real, :nan?)