From: jzakiya <jzakiya@...>
Date: 2010-01-01T08:15:05+09:00
Subject: Re: Roots Module

Added ability to now do roots Complex numbers.
After thinking about it, it was elegantly simple to do. :)

So now, module Roots can be mixed in with the numeric classes
Integer, Float, Complex, and Rational to find the roots of all
the numeric number classes, eg:

Complex(2,11).root(3) => (2.0+1.0i) [Ruby 1.9.1p243]

Today is the last day of 2009, and the penultimate day of
Kwanzaa -- Kuumba/Creativity.

May the New Year bring much needed Justice, Peace, and Love to, and
from, the Human species.

File roots.rb
--------------------------------
module Roots
=begin
Author: Jabari Zakiya,   Original: 2009-12-25
Revision-2: 2009-12-31

Mixin Roots into Integer, Float, Complex, Rational
to add methods root & roots, to find the roots of
all real and complex numbers.

Use syntax:  val.root(n,{k})
root(n,k=0) n is root 1/n exponent,  integer > 0
            k is nth ccw 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 ccw complex root    when n is even

9.root(2); -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, optional string input, are:
   0   : default (no input), return array of n ccw roots
'c'|'C': complex, return array of complex roots
'e'|'E': even, return array even numbered roots
'o'|'O': odd , return array odd  numbered roots
'i'|'I': imag, return array of imaginary roots
'r'|'R': real, return array of real roots
An empty array is returned for an opt with no members.

9334943.roots(9); -89.roots(4,'real'); 2.2.roots 3,'Im'

Can ask: How many real roots of: x.roots(n,'real').size
What's the 3rd 5th root of (4+9i): Complex(4,9).root(5,3)

Ruby currently gives 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, causing wrong root values there.
I 'fix' by clipping absolute values less than a value
I call TRIG-EPSILON so they produces the correct results.

Put trig code into separate file and 'require' into apps
to get correct|exact results for x|y axis angles and get
same accuracy for very small angles close to each 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

  require 'complex'
  include Math

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

  def roots(n,opt=0) # return array of roots, or empty
    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 = abs**n**-1
    theta = kind_of?(Complex) ? arg/n : PI/n
    roots = []
    case opt
    when /^(o|O)/  # even roots 2,4,6...
      0.step(n-1,2) {|k| roots << rootn(mag,theta,k,n)}
    when /^(e|E)/  # odd roots 1,3,5...
      1.step(n-1,2) {|k| roots << rootn(mag,theta,k,n)}
    when /^(r|R)/  # real roots     Complex(x,0) = (x+i0)
      n.times {|k|
        x=rootn(mag,theta,k,n); roots << x if x.imag == 0}
    when /^(i|I)/  # imaginry roots Complex(0,y) = (0+iy)
      n.times {|k|
        x=rootn(mag,theta,k,n); roots << x if x.real == 0}
    when /^(c|C)/  # complex roots  Complex(x,y) = (x+iy)
      n.times {|k|
        x=rootn(mag,theta,k,n); roots << x unless
        x.imag == 0 || x.real == 0}
    else           # all n roots
      n.times {|k| roots << rootn(mag,theta,k,n)}
    end
    return roots
  end

  private # don't show as methods in mixin class
  TRIG_EPSILON = 2.5e-16
  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 # better than tan

  def rootn(mag,theta,k=0,n=1) # nth root of real|complex
    angle_n = kind_of?(Complex) ? theta+(2*k*PI)/n :
              self > 0 ? 2*(k+1)*theta : (2*k+1)*theta
    mag*Complex(cosine(angle_n),sine(angle_n))
  end
end