From: jzakiya <jzakiya@...>
Date: 2009-12-26T07:15:06+09:00
Subject: Re: Roots Module

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