From: jzakiya Date: 2009-12-27T05:40:26+09:00 Subject: Re: Roots Module On Dec 26, 2:33 am, jzakiya wrote: > On Dec 26, 2:18 am, jzakiya wrote: > > > > > On Dec 25, 5:13 pm, jzakiya wrote: > > > > On Dec 25, 4:42 pm, jzakiya 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 > > >   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 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 = 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)=-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 =begin 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=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 94835805813.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 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 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 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) # 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 = 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 # 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 = 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 0 ? 2*(k+1)*theta : (2*k+1)*theta mag*Complex(cosine(angle_n),sine(angle_n)) end end