From: oinkoink+unet@... (Bret Jolly)
Date: 2004-06-10T02:33:41+09:00
Subject: Re: [RCR] Floating point division operator /. (or fdiv method)

Michael Neumann <mneumann@ntecs.de> wrote in message news:<20040605102959.GA2225@miya.intranet.ntecs.de>...
> On Sat, Jun 05, 2004 at 07:03:03PM +0900, Mark Sparshatt wrote:
> > Try this
> > 
> > >>> 1/2
>  0
> > >>> from __future__ import division
> > >>> 1/2
>  0.5
> > >>> 1//2
> > 0
> > 
> > That's using version 2.3.3
> 
> Ah thanks.... would be great to see this in Ruby too (without breaking
> existing code).

It would *not* be great to see this in Ruby: the quotient of two integers
should be a rational, not a float.  The behavior given by the mathn library
is the correct one.  If the quotient of two rationals is a float, you get
severe errors using the matrix standard library with ill-conditioned matrices.
The Hilbert matrix, which arises in the the theory of least-square fits, is
a good example.

First recall that you *must* require mathn if you use the matrix library.
Here's what happens if you forget:

irb(main):001:0> require 'matrix'
true
irb(main):002:0> m = Matrix[[2, 3], [3, 5]]
Matrix[[2, 3], [3, 5]]
irb(main):003:0> m.det
4                        # WRONG
irb(main):004:0> m.inv
Matrix[[0, 0], [-1, 0]]  # WRONG
irb(main):005:0> m*m.inv
Matrix[[-3, 0], [-5, 0]] # WRONG: m*m.inv must be the identity matrix
irb(main):006:0> require 'mathn'
true
irb(main):007:0> m.det
1                        # correct
irb(main):008:0> m.inv
Matrix[[5, -3], [-3, 2]] # correct
irb(main):009:0> m*m.inv
Matrix[[1, 0], [0, 1]]   # correct

Now let's look at what happens if you use float results for integer division.
I define two functions: hilbert and hilbert_float, which create a Hilbert
matrix and its floating-point version, respectively.
# hilbert.rb
# Hilbert matrix of order n
require 'mathn'  # ESSENTIAL!

# exact (rational entries)
def hilbert(n)
  row_array = []
  (1..n).each {|r|
     row = []
     (r..(n+r-1)).each{|s| row.push(1/s)}
     row_array.push(row)
   }
  Matrix[*row_array]
end

# floating-point entries
def hilbert_float(n)
  row_array = []
  (1..n).each {|r|
     row = []
     (r..(n+r-1)).each{|s| row.push((1/s).to_f)}
     row_array.push(row)
   }
  Matrix[*row_array]
end

class Matrix
  # find off-diagonal element of maximum magnitude
  def biggest_offdiag
    biggest = 0
    (0...row_size).each {|i|
      (0...column_size).each {|j|
        if i != j and biggest.abs < self[i, j].abs
          biggest = self[i, j]
        end
      }
    }
    biggest
  end
end
# end of file

A matrix times its inverse is the identity matrix -- 1s on the main
diagonal and 0s elsewhere.  So for any invertible matrix m,
(m*m.inv).biggest_offdiag should be 0.  Let's try it.
irb(main):001:0> require 'mathn'
true
irb(main):002:0> require 'hilbert'
true
irb(main):003:0> ((m=hilbert(15))*m.inv).biggest_offdiag
0
irb(main):004:0> ((m=hilbert_float(15))*m.inv).biggest_offdiag
-569.0       # INCREDIBLY WRONG!

The quotient of two integers is an exact rational.  To make the
quotient default to (inexact) floats can lead to appalling numerical
errors as I have just shown.

Regards, Bret