ruby-core

I have detected a bug in the implementation of Complex division.
The bug gives internal underflow or overflow which causes the result  
to contain
nans or otherwise erroneous result when the result would in fact be  
perfectly representable
with floating point real and imaginary parts.
Erroneous results will ALWAYS be returned when the absolute value of  
the denominator
exceeds the square root of the maximum representable float.
The following patch fixes this problem,and some others
Index: lib/complex.rb
===================================================================
--- lib/complex.rb	(revision 14846)
+++ lib/complex.rb	(working copy)
@@ -192,8 +192,23 @@
    #
    def / (other)
      if other.kind_of?(Complex)
-      self*other.conjugate/other.abs2
+      denominator = other.abs2
+      numerator = self * other.conjugate
+      if numerator.zero? or (denominator.kind_of?(Float) and ! 
denominator.finite?)
+        magn = other.abs
+        self * Complex(other.real / magn, - other.image / magn) / magn
+      else
+        numerator / denominator
+      end
      elsif Complex.generic?(other)
+      # If dividing with an integer, we must convert to a type that  
can divide the left operand without truncation
+      if other.kind_of?(Integer)
+         if defined?(Rational)
+           other = Rational(other,1)
+         else
+           other = Float(other)
+         end
+      end
        Complex(@real/other, @image/other)
      else
        x, y = other.coerce(self)

An incomplete unit test file follows
#!/usr/bin/ruby

#Delete this line to test with standard distribution complex.rb
$:.unshift Dir::getwd

require 'test/unit'

require 'complex'

class TestComplexBasics < Test::Unit::TestCase

   def approx_equal?(c1,c2)
     minerrlimit = 100.0 * Float::MIN
     abs1 = c1.abs
     abs2 = c2.abs
     if abs1 > abs2
       maxabs = abs1
     else
       maxabs = abs2
     end
     rel = maxabs * Float::EPSILON * 25
     if rel < minerrlimit
       rel = minerrlimit
     end
     d = (c1 - c2).abs
     return d < rel
   end

   def test_plus_integer
     assert_equal(Complex(8,7),Complex(3,4)+Complex(5,3))
     assert_equal(Complex(5,1),Complex(4,1) + 1)
     assert_equal(Complex(5,1),1+ Complex(4,1))
   end

   def test_minus_integer
   	assert_equal(Complex(2,4),Complex(3,2)-Complex(1,-2))
   	assert_equal(Complex(4,5),Complex(5,5) -1)
   	assert_equal(Complex(1,-2), 2 - Complex(1,2))
   end

   def test_times_integer
     assert_equal(Complex(-1,0),Complex(0,1) * Complex(0,1))
     assert_equal(Complex(3,3), Complex(1,1) * 3)
     assert_equal(Complex(4,-2), 2 * Complex(2,-1))
   end

   def test_plus_float
     assert(approx_equal? 
(Complex(1.0,4.0),Complex(2.0,2.0)+Complex(-1.0,2.0)))
   end

   def test_minus_float
     assert(approx_equal?(Complex(3.0,0.0),Complex(2.0,2.0)- 
Complex(-1.0,2.0)))
   end

   def test_times_float
     assert(approx_equal?(Complex(-2,1),Complex(1,2)*Complex(0,1)))
   end

   def test_division_float
     [Complex(1.0,25.0), Complex(8,97), Complex(Float::MAX/2,  
Float::MAX/3),
       Complex(100.0 / Float::MAX, 200.0 / Float::MAX)].each { | c |
       assert(approx_equal?(Complex(1.0,0.0), c / c))
     }
     assert(approx_equal?(Complex(-1.0,0.0),Complex(0.0,1.0)/ 
Complex(0.0,-1.0)))
   end

   def test_eps_log
     assert(approx_equal? 
(Complex(2.0,1.0),Math.log(Math.exp(Complex(2.0,1.0)))))
   end
end
-----------------------------------
See the amazing new SF reel: Invasion of the man eating cucumbers from  
outer space.
On congratulations for a fantastic parody, the producer replies :  
"What parody?"

Tommy Nordgren
tommy.nordgren@comhem.se

Attachments (2)

tc_complex.rb (1.78 KB, application/x-sh)

ruby-changes.patch (999 Bytes, text/x-diff)

Index: lib/complex.rb
===================================================================
--- lib/complex.rb	(revision 14846)
+++ lib/complex.rb	(working copy)
@@ -192,8 +192,23 @@
   #
   def / (other)
     if other.kind_of?(Complex)
-      self*other.conjugate/other.abs2
+      denominator = other.abs2
+      numerator = self * other.conjugate
+      if numerator.zero? or (denominator.kind_of?(Float) and !denominator.finite?)
+        magn = other.abs
+        self * Complex(other.real / magn, - other.image / magn) / magn
+      else
+        numerator / denominator
+      end
     elsif Complex.generic?(other)
+      # If dividing with an integer, we must convert to a type that can divide the left operand without truncation
+      if other.kind_of?(Integer)
+         if defined?(Rational)
+           other = Rational(other,1)
+         else
+           other = Float(other)
+         end
+      end
       Complex(@real/other, @image/other)
     else
       x, y = other.coerce(self)

Thread

Prev Next

In This Thread

Prev Next