Avoid overflow and underflow in Complex#/

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

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

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