Roots Module

As a Christmas/Holiday present to the Ruby community I finally
completed creating a versatile module which will find the accurate
real and complex roots of +|- real values.

This work evolved from concerns raised in this thread:

http://groups.google.com/group/comp.lang.ruby/browse_thread/thread/9b8a417b8f3b79f3#

from which I expanded on this work and incorporated some suggestions.

I have tested this mixed in this module to classes Float, Integer, and
BigDecimal successfully.

I tested this on Ruby versions 1.9.1p243, 1.8.7p174 and
Ruby Enterprise Edition 1.8.6 20090610.

Constructive  feedback is welcome.

File: roots.rb
------------------------------------------------------
module Roots
=begin
Mixin Roots into Integer, BigDecimal, Float to add
methods root & roots, to find the real|complex roots
of +|- real values.

Use syntax:  val.root(n,{k}}
root(n,k=0) n is root (1/n) exponent, integer > 0,
            k is nth root 1..n , integer >=0
If k not given default root returned, which are:
for +val => real root |val**(1.0/n)|
for -val => real root -|val**(1.0/n)| when n is odd
for -val => first complex root -val   when n is even
for any val => first ccw root real|complex when k > 0

9.root(2); 8.root(3,1), -32.root(5,3), -100.43.root 6,6

Use syntax:  val.roots(n,{opt}}
roots(n,opt=0) n is root (1/n) exponent, integer > 0,
               opt area optional input options below
0      : default (no input), return array of n ccw root vals
'c'|'C': complex, return array of complex roots, nil if none
'e'|'E': even, return array even numbered roots, nil if none
'o'|'O': odd , return array odd  numbered roots, nil if none
'i'|'I': imag, return array of imaginary roots, nil if none
'r'|'R': real, return array of real roots, nil if none

293481349432.roots(9); -892.roots(4,'real'); 22.2.roots 3,'i'
=end

  require 'complex'
  include Math

  def root(n,k=0) # k=1..n; k=0 for default root
    raise "Root not integer >0" unless n.kind_of?(Integer) && n>=1
	raise "Index k not an integer" unless k.kind_of?(Integer)
    return self if n == 1 || self == 0
    mag = self.abs**n**-1
	return rootn(mag,k-1,PI/n) if k>0 # nth root, k = 1..n, of any real
	return mag if self > 0            # pos real default
	return -1*mag if n&1 == 1         # neg real default, n odd
	return rootn(mag,0,PI/n)          # neg real default, n even
  end

  def roots(n,opt=0) #  returns an array of values, or nil if option
not valid
    raise "Root not integer >0" unless n.kind_of?(Integer) && n>=1
	raise "Invalid option" unless opt == 0 || opt =~ /^(c|e|i|o|r|C|E|I|O|
R)/
    return [self] if n == 1 || self == 0
    mag = self.abs**n**-1
	roots = []; theta = PI/n
	case opt
	  when /^(o|O)/  # even roots 2,4,6...
	    0.step(n-1,2) {|k| roots << rootn(mag,k,theta)}
	  when /^(e|E)/  # odd roots 1,3,5...
	    1.step(n-1,2) {|k| roots << rootn(mag,k,theta)}
	  when /^(r|R)/  # real roots  Complex(x,0) =(x+i0)
	    n.times {|k|
		  x=rootn(mag,k,theta); roots << x if x.imag == 0}
	  when /^(i|I)/  # imaginry roots  Complex(0,x) = (0+ix)
	    n.times {|k|
		  x=rootn(mag,k,theta); roots << x if x.real == 0}
	  when /^(c|C)/  # complex roots Complex(x,y) = (x+iy)
	    n.times {|k|
		  x=rootn(mag,k,theta); roots << x unless
		  x.imag == 0 || x.real == 0}
	  else
	      n.times {|k| roots << rootn(mag,k,theta)}
    end
	return roots.empty? ? nil : roots
  end

=begin
Ruby currently produces incorrect values for x|y axis angles.
cos PI/2   =>  6.12303176911189e-17
sin PI     =>  1.22460635382238e-16
cos 3*PI/2 => -1.83690953073357e-16
sin 2*PI   => -2.44921970764475e-16

These all should be 0.0, which causes incorrect root values there.
I 'fix' these errors by clipping absolute values less than a value
I call TRiG-EPSILON so they produces the correct results.

Extract this code into separate file and 'require' into your apps
to get correct|exact results for x|y axis angles and still get same
accurrcay for extremely small delta angles to the x|y axis.
cosine(89.9999999999*PI/180) =>  1.74534333112399e-11
cosine(90.0*PI/180)          =>  0.0
cosine(90.0000000001*PI/180) => -1.74543140899798e-12
=end

  protected
  TRIG_EPSILON = 1e-15
  def sine(x); y=sin(x); y.abs < TRIG_EPSILON ? 0.0:y end
  def cosine(x); y=cos(x); y.abs<TRIG_EPSILON ? 0.0:y end
  def tangent(x); sine(x)/cosine(x) end # not used here but more
correct

  def angle(k,theta) # roots 1..n --> k = 0..n-1
    angle = self > 0.0 ? 2*(k+1)*theta : (2*k+1)*theta
  end

  def rootn(mag,k,theta)
    a = angle(k,theta); mag*Complex(cosine(a),sine(a))
  end
end