金光です。

> [mailto:ruby-math-admin / ruby-lang.org] On Behalf Of 正木 功
> Sent: Tuesday, June 12, 2001 10:50 AM
> 
> 正木です。

どもっ。


> 誤差項付きの Real Class を作ってみました。

これって、どんなときどういう風に使うんですか?
ちょっとアプリを見せてもらえません?


> --------
> class Real  < Numeric
>   @@error=0
>   def initialize(seq,error=nil)
>     case error
>     when Sequence
>       @seq=seq
>       @error=error
>     when "D"
>       @seq=seq
>       @error=@seq.Diff.abs.Shift
>     when "CF"
>       @seq=CFraction_to_Sequence(seq)
>       @error=@seq.Diff.abs.Shift
>     when "GCF"
>       @seq=GCFraction_to_Sequence(seq)
>       @error=@seq.Diff.abs.Shift
>     when "S"
>       @seq=seq.Sum
>       @error=seq.abs.Shift
>     when "A"
>       @seq=seq.alt.Sum
>       @error=seq.abs.Shift
>     when "W"
>       @seq=Sequence.new([],seq){"|n,seq| seq[n][0]"}
>       @error=Sequence.new([],seq){"|n,seq| seq[n][1]"}
>     end
>     @max=@seq+@error
>     @min=@seq-@error
>   end
>   def [](n)
>     @seq[n]
>   end
>   def each
>     i=0
>     loop do
>       yield(@seq[i],@error[i],i)
>       i+=1
>     end
>   end
>   def find
>     r=nil;each{|x,e,i| if yield(x,e,i) then r=[x,e,i];break;end};r
>   end
>   def +(a)
>     case a
>     when Real,Integer,Rational,Float
>       y=Real(a)
>       Real.new(@seq+y.seq,@error+y.error)
>     else
>       x , y = a.coerce(self)
>       x + y
>     end
>   end
>   def -(a)
>     case a
>     when Real,Integer,Rational,Float
>       y=Real(a)
>       Real.new(@seq-y.seq,@error+y.error)
>     else
>       x , y = a.coerce(self)
>       x - y
>     end
>   end
>   def	-@
>     Real.new(-@seq,@error)
>   end
>   def abs
>     Real.new(@seq.abs,@error)
>   end
>   def absmax
>     @seq.abs+@error
>   end
>   def absmin
>     x=self
>     Sequence.new([],x){"|n,x|  [x[n].abs-x.error[n],0].max "}
>   end
>   def *(y)
>     case y
>     when Integer,Rational
>       Real.new(@seq*y,@error*y.abs)
>     when Float
>       self*Real(y)
>     when Real
>       Real.new(@seq*y.seq,absmax*y.absmax- / seq.abs*y.seq.abs)
>     else
>       x , y = y.coerce(self)
>       x*y
>     end
>   end
>   def /(y)
>     case y
>     when Integer,Rational
>       Real.new(@seq/y,@error/y.abs)
>     when Float
>       self/Real(y)
>     when Real
>       Real.new(@seq/y.seq,absmax/y.absmin- / seq.abs/y.seq.abs)
>     else
>       x , y = y.coerce(self)
>       x/y
>     end
>   end
>   def **(m)
>     Real.new(@seq**m,absmax**m- / seq.abs**m)
>   end
>   def inverse
>     error=absmin.inverse- / seq.abs.inverse
>     Real.new(@seq.inverse,error)
>   end
>   def approx(eps=@@error)
>     find{|x,e,i| e <= eps }[0]
>   end
>   def Shift(m=1)
>     Real.new(@seq.Shift(m),@error.Shift(m))
>   end
>   def printd(n=16,r=10)
>     approx(r**(-n)).printd(n,r)
>   end
>   def <=> (other)
>     case other
>     when Integer,Rational,Float,Real
>       x=self-Real(other)
>       n=0
>       while true
> 	return 1 if x.min[n]>0
> 	return -1 if x.max[n]<0
> 	return 0 if x.absmax[n] <= @@error
> 	n+=1
>       end
>     else
>       x , y = other.coerce(self)
>       return x <=> y
>     end
>   end
>   def coerce(other)
>     return Real(other),self
>   end
>   def Real.error=(e)
>     @@error=e
>   end
>   def Real.error
>     @@error
>   end
>   def sqrt
>     x=self
>     s=Sequence.new([Real.sqrt(x[0])[0]],x){"|n,x|
> Real.sqrt(x[n],self[n-1])[0]
> "}    error=Sequence.new([],x,s){"|n,x,s| 
> [Real.sqrt(x.max[n],s[n]).max[0]-s[n
> ],s[n]-Real.sqrt(x.absmin[n],s[n]).min[0]].max"}
>     Real.new(s,error)
>   end
>   def exp
>     x=self
>     s=Sequence.new([],x){"|n,x| Real.exp(x[n])[n]"}
>     e=Sequence.new([],x,s){"|n,x,s| Real.exp(x.max[n]).max[n]-s[n]"}
>     Real.new(s,e)
>   end
>   def log
>     x=self
>     s=Sequence.new([],x){"|n,x| Real.log(x[n])[n]"}
>     e=Sequence.new([],x,s){"|n,x,s| Real.log(x.max[n]).max[n]-s[n]"}
>     Real.new(s,e)
>   end
>   def cos
>     pi=Realpi
>     @@error=1/8
>     x=self
>     return (-x).cos if x<0
>     return -(x-pi).cos if x>pi
>     return -(x-pi/2).sin if x>pi/2
>     return (pi/2-x).sin if x>pi/4
>     s=Sequence.new([],x){"|n,x| Real.cos(x[n])[n]"}
>     e=Sequence.new([],x,s){"|n,x,s| 
> [Real.cos(x.absmin[n]).max[n]-s[n],s[n]-Re
> al.cos(x.absmax[n]).min[n]].max if x.min[n]>-2 && x.max[n]<2 "}
>     Real.new(s,e)
>   end
>   def sin
>     pi=Realpi
>     @@error=1/8
>     x=self
>     return -(-x).sin if x<0
>     return -(x-pi).sin if x>pi
>     return (x-pi/2).cos if x>pi/2
>     return (pi/2-x).cos if x>pi/4
>     s=Sequence.new([],x){"|n,x| Real.sin(x[n])[n]"}
>     e=Sequence.new([],x,s){"|n,x,s| 
> [Real.sin(x.max[n]).max[n]-s[n],s[n]-Real.
> sin(x.min[n]).min[n]].max if x.max[n] < 3/2 && x.min[n] > -3/2"}
>     Real.new(s,e)
>   end
>   def tan
>     sin/cos
>   end
>   attr_reader :seq,:error,:max,:min
> end
> class Object
>   def Real(x)
>     case x
>     when Real
>       x
>     when Integer,Rational
>       Real.new(Sequence.new([],x){"|n,x| x"},Sequence.new(){"|n| 0"})
>     when Float
>       x.to_Real
>     end
>   end
> end
> class Fixnum,Bignum
>   alias / rdiv
> end
> def GCFraction_to_Sequence(s)
>   a=Sequence.new([0,1],s){"|k,s| p,q=s[k-2];p*self[k-2]+q*self[k-1] "}
>   b=Sequence.new([1,0],s){"|k,s| p,q=s[k-2];p*self[k-2]+q*self[k-1] "}
>   (a/b).Shift(2)
> end
> def Real.sqrt(x,a=1)
>   s=Sequence.new([(a+x/a)/2],x){"|n,x| a=self[n-1];(a+x/a)/2 "}
>   error=Sequence.new([],x,s){"|n,x,s| y=s[n];y-x/y "}
>   Real.new(s,error)
> end
> def Real.exp(x)
>   if x<0
>     Real.new(ExpS[x],"S")
>   elsif x>Rational(1,2)
>     Real.exp(x/2)**2
>   else
>     Real.exp(-x).inverse
>   end
> end
> def Real.e
>   Real.exp(1)
> end
> def Real.cos(x)
>   Real.new(Sequence.new([],x){"|n,x| ExpS[x][2*n]"},"A")
> end
> def Real.sin(x)
>   Real.new(Sequence.new([],x){"|n,x| ExpS[x][2*n+1]"},"A")
> end
> def Real.atan(x)
>   return -Real.atan(-x) if x<0
>   return Realpi/2-Real.atan(1/x) if x>1
>   Real.new(OddPower[x]/Odd,"A")
> end
> def Real.pi
>   16*Real.atan(1/5)-4*Real.atan(1/239)
> end
> Realpi=Real.pi
> def Real.log(x)
>   return -(Real.log(1/x)) if x<1
>   y=(x-1)/(x+1)
>   term=OddPower[y]/Odd*2
>   Real.new(term.Sum,term.Shift*(x+3)/4)
> end
> 
> --------
> 
> 問題点がいくつかあります。
> 1、本当に有理数の Cauchy 列かどうかの判定をしない。
> 2、誤差項が正しいかどうかの判定もしない。
> 3、収束の保証されている手続きによって得られる実数のみを扱うことにして
> も、
> それは実数全体から見ればごく一部に過ぎない。従って class 名が不適当。
> (computer で扱える数は可算個しかない筈なので、本当の実数全体の class
が
> 作れると思う人もいないでしょうが)
> 4、2つの実数が本当に等しいかどうかの判定はできない。
> 2つの実数が等しいときに @@error=0 として <=> を実行すると無限 loop
> に陥ってしまう。
> 等しくない事が分かっている場合、大小比較は可能。例えば
> p Math::PI < Real.pi #=> true
> 
> 

M.Kanemitsu
―――――――――――――――
金光雅夫 masao-k / a-net.email.ne.jp

http://www.ne.jp/asahi/masao-k/home/
http://www14.freeweb.ne.jp/art/soshikon/
http://www15.freeweb.ne.jp/computer/ruby256/
〒216-0031 川崎市宮前区神木本町5-14-12
自宅: 044-877-5006
携帯: 090-2753-5292