################################################################## ### # # File: Apc.rb # # Subject: Class for arbitrary precision complex. # # Author: Dennis J. Darland # # Date: April 2, 2007 # ############################################################################ # # Copyright (C) 2007 Dennis J. Darland# # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # ############################################################################ class Apc # # @re is real comonent # @im is imaginary component def initialize(re,im) @re = re @im = im end def re @re end def im @im end def setcconst(c) @@cconst = c end def setrconst(c) @@rconst = c end def neg Apc.new(@re.neg,@im.neg) end def to_s if @im > @@rconst.zero then "#{@re.to_s} + #{@im.to_s}I" else "#{@re.to_s} - #{@im.neg.to_s}I" end end def display_val(label) if RAILS_MODE == 0 then puts label + " = " + self.to_s end end def +(other) Apc.new(@re+other.re,@im+other.im) end def -(other) Apc.new(@re-other.re,@im-other.im) end def *(other) Apc.new(@re*other.re-@im*other.im,@im*other.re+@re*other.im) end def /(other) de = other.re * other.re + other.im * other.im Apc.new((@re * other.re + @im * other.im)/de ,(@im * other.re - @re * other.im)/de) end def ==(other) (@re == other.re) && (@im == other.im) end def eq(other) (@re.eq(other.re)) && (@im.eq(other.im)) end def abs if self == @@cconst.zero return @@rconst.zero end (@re * @re + @im * @im).sqrt end def exp w = @re.exp Apc.new(w * @im.cos,w * @im.sin) end def sin y1 = (self * @@cconst.i).exp y2 = (self * @@cconst.minus_i).exp (y1-y2) / @@cconst.two_i end def cos y1 = (self * @@cconst.i).exp y2 = (self * @@cconst.minus_i).exp (y1+y2) / @@cconst.two end def tan self.sin/self.cos end def sinh y1 = self.exp y2 = (self * @@cconst.minus_one).exp (y1-y2) / @@cconst.two end def cosh y1 = self.exp y2 = (self * @@cconst.minus_one).exp (y1+y2) / @@cconst.two end def tanh self.sinh/self.cosh end def log dd = self.abs # puts "cx log dd = " + dd.to_s if dd > @@rconst.zero then lg = dd.log # puts "cx log lg = dd.log = " + lg.to_s else lg = @@rconst.zero #really error! end if self.re == @@rconst.zero then if self.im > @@rconst.zero then ang = @@rconst.pi/@@rconst.two else ang = @@rconst.pi/@@rconst.minus_two end else ang = self.calc_lg_ang end return Apc.new(lg,ang) end def calc_lg_ang if @re >= @@rconst.zero && @im >= @@rconst.zero then return (@im / @re).atan elsif @re >= @@rconst.zero && @im < @@rconst.zero then return (@im / @re).atan elsif @re <= @@rconst.zero && @im >= @@rconst.zero then return @@rconst.pi + (@im / @re).atan else return @@rconst.pi * @@rconst.minus_one + (@im / @re).atan end end def log10 self.log/@@cconst.ten.log end def **(other) a = self.clone unless other.kind_of?(Apc) #other is integer b = @@cconst.one if other > 0 then other.times do b = b.clone * a end elsif other < 0 then (-other).times do b = b.clone * a end b = @@cconst.one/b.clone end return b else if other == @@cconst.zero then return @@cconst.one elsif other == @@cconst.one then return a elsif a == @@cconst.zero then return a elsif a == @@cconst.one then return a else lg = a.log lgt = lg * other b = lgt.exp return b end end end def sqrt self ** @@cconst.half end # def **(other) # a = self.clone # unless other.kind_of?(Apc) # b = @@cconst.one.clone # if other > 0 then # other.times do # b = a * b.clone # end # return b # elsif other < 0 # (-other).times do # b = a * b.clone # end # c = @@cconst.one.clone/b.clone # return c # end # return b # else # b = other # if b == (@@cconst.zero) then return @@cconst.one end # if b == @@cconst.one then return a end # if a == @@cconst.zero then return @@cconst.zero end # return (a.log*b).exp # end # end def asin if self.abs < c(0.000001) then return self end @@cconst.minus_i * (@@cconst.i * self+(@@cconst.one - self*self).sqrt).log end def acos Apc.new(@@rconst.pi / @@rconst.two,@@rconst.zero) - self.asin end def atan if self.abs < c(0.000001) then return self end (@@cconst.i / @@cconst.two) * ((@@cconst.one - self * @@cconst.i).log - (@@cconst.one + self * @@cconst.i).log) end def asinh (self + (self * self + @@cconst.one).sqrt).log end def acosh (self + (self * self - @@cconst.one).sqrt).log end def atanh @@cconst.half * ((@@cconst.one + self).log - (@@cconst.one - self).log) end end #class Apc class ApcConst def initialize(rConst) @ZERO = Apc.new(rConst.zero,rConst.zero) @ONE = Apc.new(rConst.one,rConst.zero) @TWO = Apc.new(rConst.two,rConst.zero) @TEN = Apc.new(rConst.ten,rConst.zero) @HALF = Apc.new(rConst.half,rConst.zero) @I = Apc.new(rConst.zero,rConst.one) @MINUS_I = Apc.new(rConst.zero,rConst.minus_one) @MINUS_ONE = Apc.new(rConst.minus_one,rConst.zero) @TWO_I = Apc.new(rConst.zero,rConst.two) end #initialize def zero @ZERO end def one @ONE end def two @TWO end def i @I end def minus_i @MINUS_I end def minus_one @MINUS_ONE end def two_i @TWO_I end def ten @TEN end def half @HALF end end #ApcConst