Re: Roots Module

On Dec 26, 2:33=A0am, jzakiya <[e-mail]> wrote:
> On Dec 26, 2:18=A0am, jzakiya <[e-mail]> wrote:
>
>
>
> > On Dec 25, 5:13=A0pm, jzakiya <[e-mail]> wrote:
>
> > > On Dec 25, 4:42=A0pm, jzakiya <[e-mail]> 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/9=
b...
>
> > > > from which I expanded on this work and incorporated some suggestion=
s.
>
> > > > 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 =A0feedback is welcome.
>
> > > > File: roots.rb
> > > > ------------------------------------------------------
> > > > module Roots
> > > > =3Dbegin
> > > > Mixin Roots into Integer, BigDecimal, Float to add
> > > > methods root & roots, to find the real|complex roots
> > > > of +|- real values.
>
> > > > Use syntax: =A0val.root(n,{k}}
> > > > root(n,k=3D0) n is root (1/n) exponent, integer > 0,
> > > > =A0 =A0 =A0 =A0 =A0 =A0 k is nth root 1..n , integer >=3D0
> > > > If k not given default root returned, which are:
> > > > for +val =3D> real root |val**(1.0/n)|
> > > > for -val =3D> real root -|val**(1.0/n)| when n is odd
> > > > for -val =3D> first complex root -val =A0 when n is even
> > > > for any val =3D> 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: =A0val.roots(n,{opt}}
> > > > roots(n,opt=3D0) n is root (1/n) exponent, integer > 0,
> > > > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0opt area optional input options belo=
w
> > > > 0 =A0 =A0 =A0: 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 =A0numbered 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'
> > > > =3Dend
>
> > > > =A0 require 'complex'
> > > > =A0 include Math
>
> > > > =A0 def root(n,k=3D0) # k=3D1..n; k=3D0 for default root
> > > > =A0 =A0 raise "Root not integer >0" unless n.kind_of?(Integer) && n=
>=3D1
> > > > =A0 =A0 =A0 =A0 raise "Index k not an integer" unless k.kind_of?(In=
teger)
> > > > =A0 =A0 return self if n =3D=3D 1 || self =3D=3D 0
> > > > =A0 =A0 mag =3D self.abs**n**-1
> > > > =A0 =A0 =A0 =A0 return rootn(mag,k-1,PI/n) if k>0 # nth root, k =3D=
 1..n, of any real
> > > > =A0 =A0 =A0 =A0 return mag if self > 0 =A0 =A0 =A0 =A0 =A0 =A0# pos=
 real default
> > > > =A0 =A0 =A0 =A0 return -1*mag if n&1 =3D=3D 1 =A0 =A0 =A0 =A0 # neg=
 real default, n odd
> > > > =A0 =A0 =A0 =A0 return rootn(mag,0,PI/n) =A0 =A0 =A0 =A0 =A0# neg r=
eal default, n even
> > > > =A0 end
>
> > > > =A0 def roots(n,opt=3D0) # =A0returns an array of values, or nil if=
 option
> > > > not valid
> > > > =A0 =A0 raise "Root not integer >0" unless n.kind_of?(Integer) && n=
>=3D1
> > > > =A0 =A0 =A0 =A0 raise "Invalid option" unless opt =3D=3D 0 || opt =
=3D~ /^(c|e|i|o|r|C|E|I|O|
> > > > R)/
> > > > =A0 =A0 return [self] if n =3D=3D 1 || self =3D=3D 0
> > > > =A0 =A0 mag =3D self.abs**n**-1
> > > > =A0 =A0 =A0 =A0 roots =3D []; theta =3D PI/n
> > > > =A0 =A0 =A0 =A0 case opt
> > > > =A0 =A0 =A0 =A0 =A0 when /^(o|O)/ =A0# even roots 2,4,6...
> > > > =A0 =A0 =A0 =A0 =A0 =A0 0.step(n-1,2) {|k| roots << rootn(mag,k,the=
ta)}
> > > > =A0 =A0 =A0 =A0 =A0 when /^(e|E)/ =A0# odd roots 1,3,5...
> > > > =A0 =A0 =A0 =A0 =A0 =A0 1.step(n-1,2) {|k| roots << rootn(mag,k,the=
ta)}
> > > > =A0 =A0 =A0 =A0 =A0 when /^(r|R)/ =A0# real roots =A0Complex(x,0) =
=3D(x+i0)
> > > > =A0 =A0 =A0 =A0 =A0 =A0 n.times {|k|
> > > > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 x=3Drootn(mag,k,theta); roots <=
< x if x.imag =3D=3D 0}
> > > > =A0 =A0 =A0 =A0 =A0 when /^(i|I)/ =A0# imaginry roots =A0Complex(0,=
x) =3D (0+ix)
> > > > =A0 =A0 =A0 =A0 =A0 =A0 n.times {|k|
> > > > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 x=3Drootn(mag,k,theta); roots <=
< x if x.real =3D=3D 0}
> > > > =A0 =A0 =A0 =A0 =A0 when /^(c|C)/ =A0# complex roots Complex(x,y) =
=3D (x+iy)
> > > > =A0 =A0 =A0 =A0 =A0 =A0 n.times {|k|
> > > > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 x=3Drootn(mag,k,theta); roots <=
< x unless
> > > > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 x.imag =3D=3D 0 || x.real =3D=
=3D 0}
> > > > =A0 =A0 =A0 =A0 =A0 else
> > > > =A0 =A0 =A0 =A0 =A0 =A0 =A0 n.times {|k| roots << rootn(mag,k,theta=
)}
> > > > =A0 =A0 end
> > > > =A0 =A0 =A0 =A0 return roots.empty? ? nil : roots
> > > > =A0 end
>
> > > > =3Dbegin
> > > > Ruby currently produces incorrect values for x|y axis angles.
> > > > cos PI/2 =A0 =3D> =A06.12303176911189e-17
> > > > sin PI =A0 =A0 =3D> =A01.22460635382238e-16
> > > > cos 3*PI/2 =3D> -1.83690953073357e-16
> > > > sin 2*PI =A0 =3D> -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) =3D> =A01.74534333112399e-11
> > > > cosine(90.0*PI/180) =A0 =A0 =A0 =A0 =A0=3D> =A00.0
> > > > cosine(90.0000000001*PI/180) =3D> -1.74543140899798e-12
> > > > =3Dend
>
> > > > =A0 protected
> > > > =A0 TRIG_EPSILON =3D 1e-15
> > > > =A0 def sine(x); y=3Dsin(x); y.abs < TRIG_EPSILON ? 0.0:y end
> > > > =A0 def cosine(x); y=3Dcos(x); y.abs<TRIG_EPSILON ? 0.0:y end
> > > > =A0 def tangent(x); sine(x)/cosine(x) end # not used here but more
> > > > correct
>
> > > > =A0 def angle(k,theta) # roots 1..n --> k =3D 0..n-1
> > > > =A0 =A0 angle =3D self > 0.0 ? 2*(k+1)*theta : (2*k+1)*theta
> > > > =A0 end
>
> > > > =A0 def rootn(mag,k,theta)
> > > > =A0 =A0 a =3D angle(k,theta); mag*Complex(cosine(a),sine(a))
> > > > =A0 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:
>
> > =A0 def angle_n(k,theta) # roots 1..n --> k =3D 0..n-1
> > =A0 =A0 self > 0 ? 2*(k+1)*theta : (2*k+1)*theta
> > =A0 end
>
> > =A0 def rootn(mag,k,theta) # nth root of a real value
> > =A0 =A0 a =3D angle_n(k,theta); mag*Complex(cosine(a),sine(a))
> > =A0 end
>
> Actually, just eliminate 'def angle_n' since it was only
> use in rootn, and just have rootn be this:
>
> =A0 def rootn(mag,k,theta) # nth root of a real value
> =A0 =A0 angle_n =3D self > 0 ? 2*(k+1)*theta : (2*k+1)*theta
> =A0 =A0 =A0 =A0 mag*Complex(cosine(angle_n),sine(angle_n))
> =A0 end

OK, here's the final (for now) Roots module, all
cleaned up and tweaked.

New changes:
I made roots method return just an array in all cases.
So now, if an option doesn't have members, it returns
an empty array, [], instead of nil.

This allows you to now consistently do things like:
How many real roots for x? -- x.roots(n,'real').size
If there are none, you get '0', instead of exception
error of NoMethod for Nilclass for size.

I also set TRIG-EPSILON =3D 2.5e-16 to get it as close
to the largest incorrect value for cos|sin for the x|y
axis angles, which is sin(2*PI)=3D-2.44921970764475e-16
on my Intel P4 cpu system. Set for equivalent results
for your cpu (32|64-bit) environment (AMD, PPC, etc).

I also noticed when Roots is mixedin to Integer it will
also take care of Bignums too, but you have to mixin to
Float separately to take care of them too.

Finally, for those who don't know|care|believe this
stuff matters, or has any uses, check out:

Project Euler -- http://projecteuler.net/

Ruby users and aficionados have a pretty strong throng
there (along with other languages), and to get the
"correct answers" exactness matters.

So here is the current final Roots module, cleaned up
with (hopefully) no typos, and formatted here to cleanly
fit within the column limitations for posting here.

If/when I get the desire, I may create a gem out of it,
with DOC markup, and all the rest.

Finally, today is December 26, the first day of Kwanzaa,
which celebrates the principle of Umoja -- unity.

Peace and Blessings to All.

File roots.rb
--------------------------------
module Roots
=3Dbegin
Author: Jabari Zakiya,   Original: 2009-12-25
Revision-1: 2009-12-26

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=3D0) n is root 1/n exponent,  integer > 0
            k is nth ccw root 1..n , integer >=3D0
If k not given default root returned, which are:
for +val =3D> real root  |val**(1.0/n)|
for -val =3D> real root -|val**(1.0/n)| when n is odd
for -val =3D> first ccw complex root    when n is even

94835805813.root(2); -32.root(5,3), -100.43.root 6,6

Use syntax:  val.roots(n,{opt})
roots(n,opt=3D0) 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

Ruby currently gives incorrect values for x|y axis angles
cos PI/2   =3D>  6.12303176911189e-17
sin PI     =3D>  1.22460635382238e-16
cos 3*PI/2 =3D> -1.83690953073357e-16
sin 2*PI   =3D> -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) =3D>  1.74534333112399e-11
cosine(90.0*PI/180)          =3D>  0.0
cosine(90.0000000001*PI/180) =3D> -1.74543140899798e-12
=3Dend

  require 'complex'
  include Math

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

  def roots(n,opt=3D0) # return array of roots, or empty
    raise "Root not integer >0" unless
                           n.kind_of?(Integer) && n>=3D1
    raise "Invalid option" unless opt =3D=3D 0 ||
                           opt =3D~ ^(c|e|i|o|r|C|E|I|O|R)/
    return [self] if n =3D=3D 1 || self =3D=3D 0
    mag =3D self.abs**n**-1
    roots =3D []; theta =3D 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) =3D(x+i0)
      n.times {|k|
        x=3Drootn(mag,k,theta); roots << x if x.imag =3D=3D 0}
    when /^(i|I)/  # imaginry roots Complex(0,x) =3D (0+ix)
      n.times {|k|
        x=3Drootn(mag,k,theta); roots << x if x.real =3D=3D 0}
    when /^(c|C)/  # complex roots  Complex(x,y) =3D (x+iy)
      n.times {|k|
        x=3Drootn(mag,k,theta); roots << x unless
        x.imag =3D=3D 0 || x.real =3D=3D 0}
    else           # all n roots
      n.times {|k| roots << rootn(mag,k,theta)}
    end
    return roots
  end

  private # don't show as methods in mixin class
  TRIG_EPSILON =3D 2.5e-16
  def sine(x); y=3Dsin(x); y.abs < TRIG_EPSILON ? 0.0:y end
  def cosine(x); y=3Dcos(x); y.abs<TRIG_EPSILON ? 0.0:y end
  def tangent(x); sine(x)/cosine(x) end # better than tan

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