Avoid overflow and underflow in Complex#/
MODIFY
mrbgems/mruby-complex/mrblib/complex.rb
mrbgems/mruby-complex/mrblib/complex.rb
@@ -45,8 +45,7 @@ class Complex < Numeric
def /(rhs)
if rhs.is_a? Complex
- div = rhs.real * rhs.real + rhs.imaginary * rhs.imaginary
- Complex((real * rhs.real + imaginary * rhs.imaginary) / div, (rhs.real * imaginary - real * rhs.imaginary) / div)
+ __div__(rhs)
elsif rhs.is_a? Numeric
Complex(real / rhs, imaginary / rhs)
end

MODIFY
mrbgems/mruby-complex/src/complex.c
mrbgems/mruby-complex/src/complex.c
#include <mruby.h>
#include <mruby/class.h>
#include <mruby/numeric.h>
+#include <math.h>
#ifdef MRB_WITHOUT_FLOAT
# error Complex conflicts 'MRB_WITHOUT_FLOAT' configuration in your 'build_config.rb'
@@ -11,6 +12,12 @@ struct mrb_complex {
mrb_float imaginary;
};
+#ifdef MRB_USE_FLOAT
+#define F(x) x##f
+#else
+#define F(x) x
+#endif
+
#if defined(MRB_64BIT) || defined(MRB_USE_FLOAT)
#define COMPLEX_USE_ISTRUCT
@@ -124,6 +131,93 @@ complex_to_c(mrb_state *mrb, mrb_value self)
return self;
}
+/* Arithmetic on (significand, exponent) pairs avoids premature overflow in
+ complex division */
+struct float_pair {
+ mrb_float s;
+ int x;
+};
+
+static void
+add_pair(struct float_pair *s, struct float_pair const *a,
+ struct float_pair const *b)
+{
+ if (b->s == 0.0F) {
+ *s = *a;
+ } else if (a->s == 0.0F) {
+ *s = *b;
+ } else if (a->x >= b->x) {
+ s->s = a->s + F(ldexp)(b->s, b->x - a->x);
+ s->x = a->x;
+ } else {
+ s->s = F(ldexp)(a->s, a->x - b->x) + b->s;
+ s->x = b->x;
+ }
+}
+
+static void
+mul_pair(struct float_pair *p, struct float_pair const *a,
+ struct float_pair const *b)
+{
+ p->s = a->s * b->s;
+ p->x = a->x + b->x;
+}
+
+static void
+div_pair(struct float_pair *q, struct float_pair const *a,
+ struct float_pair const *b)
+{
+ q->s = a->s / b->s;
+ q->x = a->x - b->x;
+}
+
+static mrb_value
+complex_div(mrb_state *mrb, mrb_value self)
+{
+ mrb_value rhs;
+ struct mrb_complex *a, *b;
+ struct float_pair ar, ai, br, bi;
+ struct float_pair br2, bi2;
+ struct float_pair div;
+ struct float_pair ar_br, ai_bi;
+ struct float_pair ai_br, ar_bi;
+ struct float_pair zr, zi;
+
+ mrb_get_args(mrb, "o", &rhs);
+ a = complex_ptr(mrb, self);
+ b = complex_ptr(mrb, rhs);
+
+ /* Split floating point components into significand and exponent */
+ ar.s = F(frexp)(a->real, &ar.x);
+ ai.s = F(frexp)(a->imaginary, &ai.x);
+ br.s = F(frexp)(b->real, &br.x);
+ bi.s = F(frexp)(b->imaginary, &bi.x);
+
+ /* Perform arithmetic on (significand, exponent) pairs to produce
+ the result: */
+
+ /* the divisor */
+ mul_pair(&br2, &br, &br);
+ mul_pair(&bi2, &bi, &bi);
+ add_pair(&div, &br2, &bi2);
+
+ /* real component */
+ mul_pair(&ar_br, &ar, &br);
+ mul_pair(&ai_bi, &ai, &bi);
+ add_pair(&zr, &ar_br, &ai_bi);
+ div_pair(&zr, &zr, &div);
+
+ /* imaginary component */
+ mul_pair(&ai_br, &ai, &br);
+ mul_pair(&ar_bi, &ar, &bi);
+ ar_bi.s = -ar_bi.s;
+ add_pair(&zi, &ai_br, &ar_bi);
+ div_pair(&zi, &zi, &div);
+
+ /* assemble the result */
+ return complex_new(mrb, F(ldexp)(zr.s, zr.x), F(ldexp)(zi.s, zi.x));
+}
+
void mrb_mruby_complex_gem_init(mrb_state *mrb)
{
struct RClass *comp;
@@ -146,6 +240,7 @@ void mrb_mruby_complex_gem_init(mrb_state *mrb)
mrb_define_method(mrb, comp, "to_f", complex_to_f, MRB_ARGS_NONE());
mrb_define_method(mrb, comp, "to_i", complex_to_i, MRB_ARGS_NONE());
mrb_define_method(mrb, comp, "to_c", complex_to_c, MRB_ARGS_NONE());
+ mrb_define_method(mrb, comp, "__div__", complex_div, MRB_ARGS_REQ(1));
}
void

MODIFY
mrbgems/mruby-complex/test/complex.rb
mrbgems/mruby-complex/test/complex.rb
@@ -59,6 +59,15 @@ assert 'Complex#/' do
assert_complex Complex(-2, 9) / Complex(-9, 2), ((36 / 85) - (77i / 85))
assert_complex Complex(9, 8) / 4 , ((9 / 4) + 2i)
assert_complex Complex(20, 9) / 9.8 , (2.0408163265306123 + 0.9183673469387754i)
+ if 1e39.infinite? then
+ # MRB_USE_FLOAT in effect
+ ten = 1e21
+ one = 1e20
+ else
+ ten = 1e201
+ one = 1e200
+ end
+ assert_complex Complex(ten, ten) / Complex(one, one), Complex(10.0, 0.0)
end
assert 'Complex#==' do
@@ -76,8 +85,8 @@ assert 'Complex#abs' do
else
exp = 1021
end
- assert_true Complex(3.0*2**exp, 4.0*2**exp).abs.finite?
- assert_float Complex(3.0*2**exp, 4.0*2**exp).abs, 5.0*2**exp
+ assert_true Complex(3.0*2.0**exp, 4.0*2.0**exp).abs.finite?
+ assert_float Complex(3.0*2.0**exp, 4.0*2.0**exp).abs, 5.0*2.0**exp
end
assert 'Complex#abs2' do

Avoid overflow and underflow in Complex#/