From: jzakiya <jzakiya@...>
Date: 2009-12-26T16:35:05+09:00
Subject: Re: Roots Module

On Dec 26, 2:18 am, jzakiya <jzak...@gmail.com> wrote:
> On Dec 25, 5:13 pm, jzakiya <jzak...@mail.com> wrote:
>
>
>
> > On Dec 25, 4:42 pm, jzakiya <jzak...@mail.com> wrote:
>
> > > 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/9b...
>
> > > 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
>
> > Corrections:
> > change access control from 'protected' to 'private'
> > as those definitions shouldn't show up as methods in
> > the mixedin classes.
>
> > Also put at top of comment:
> > Author: Jabari Zakiya, 2009-12-25
>
> Correction for method conflict in Float:
>
> In 1.9.1 there is conflict in Float with method 'angle'
> thus, rename my private definitions as such:
>
>   def angle_n(k,theta) # roots 1..n --> k = 0..n-1
>     self > 0 ? 2*(k+1)*theta : (2*k+1)*theta
>   end
>
>   def rootn(mag,k,theta) # nth root of a real value
>     a = angle_n(k,theta); mag*Complex(cosine(a),sine(a))
>   end

Actually, just eliminate 'def angle_n' since it was only
use in rootn, and just have rootn be this:

  def rootn(mag,k,theta) # nth root of a real value
    angle_n = self > 0 ? 2*(k+1)*theta : (2*k+1)*theta
	mag*Complex(cosine(angle_n),sine(angle_n))
  end