From ruby-math-admin@ruby-lang.org  Fri Aug 31 17:43:56 2001
Received: from tonton.nagaokaut.ac.jp (tonton.nagaokaut.ac.jp [133.44.2.115])
	by blade.nagaokaut.ac.jp (8.8.8/8.8.8/Debian/GNU) with ESMTP id RAA20799
	for <poffice@blade.nagaokaut.ac.jp>; Fri, 31 Aug 2001 17:43:55 +0900
Resent-From: ruby-math-admin@ruby-lang.org
Received: (from root@localhost)
	by tonton.nagaokaut.ac.jp (8.11.3/8.11.3) id f7V8b6a51281
	for poffice@blade.nagaokaut.ac.jp; Fri, 31 Aug 2001 17:37:06 +0900 (JST)
	(envelope-from ruby-math-admin@ruby-lang.org)
Received: from voscc.nagaokaut.ac.jp (voscc.nagaokaut.ac.jp [133.44.1.100])
	by tonton.nagaokaut.ac.jp (8.11.3/8.11.3av) with ESMTP id f7V8b5E51274
	for <poffice@blade.nagaokaut.ac.jp>; Fri, 31 Aug 2001 17:37:05 +0900 (JST)
	(envelope-from ruby-math-admin@ruby-lang.org)
Received: from helium.ruby-lang.org (helium.ruby-lang.org [210.251.121.214])
	by voscc.nagaokaut.ac.jp (8.9.3/3.7W) id RAA59681;
	Fri, 31 Aug 2001 17:37:11 +0900 (JST)
Received: from hoyogw.ruby-lang.org (localhost [127.0.0.1])
	by helium.ruby-lang.org (Postfix) with ESMTP id A34FB393A
	for <poffice@blade.nagaokaut.ac.jp>; Fri, 31 Aug 2001 17:36:40 +0900 (JST)
Date: Wed, 29 Aug 2001 07:31:09 +0900
From: =?ISO-2022-JP?B?GyRCQDVMWiEhOHkbKEI=?= <GEC01122@nifty.ne.jp>
Reply-To: ruby-math@ruby-lang.org
Subject: [ruby-math:00612] Re: class Real
To: ruby-math@ruby-lang.org
Message-Id: <200108282231.HAA25869@ums507.nifty.ne.jp>
X-ML-Name: ruby-math
X-Mail-Count: 00612
X-MLServer: fml [fml 3.0pl#17]; post only (only members can post)
X-ML-Info: If you have a question, send e-mail with the body
	"help" (without quotes) to the address ruby-math-ctl@ruby-lang.org;
	help=<mailto:ruby-math-ctl@ruby-lang.org?body=help>
Mime-Version: 1.0
Content-Type: text/plain; charset="iso-2022-jp"
Precedence: bulk
Resent-To: poffice@blade.nagaokaut.ac.jp
Resent-Date: Fri, 31 Aug 2001 17:36:40 +0900
Resent-Message-Id: <200108311736.FMLAAA22669.ruby-math@ruby-lang.org>
X-Virus-Scanned: by AMaViS perl-10

正木です。

[ruby-math:00521] の Atan 等を以下のように変更します。
Atan は Taylor 級数に ε-algorithm を適用したものが非常に早いので
次のものに換えることにします。

Atan=Function.new({0=>0}){"|x| 
      (x<0)? -Atan[-x]:  
      (x>1)? 2*Atan[1]-Atan[1/x] : 
       RealComplex.new((OddPower[x]/Odd).alt.Sum,\"E\")"}

Log も次のものに変更します。

Log=Function.new({1=>0}){"|z| 
  if z.real>=0
    RealComplex.new((2*OddPower[(z-1)/(z+1)]/Odd).Sum,\"E\")
  else
    if z.imag>=0
      Log[-z]+Log[Complex::I]*2
    else 
      Log[-z]-Log[Complex::I]*2
    end
  end
  "}

Exp も誤差項が保証できない最初の数項を捨てる次のものが最も早いようです。

Exp=Function.new({0=>1}){"|x| n=Integer(2*x.abs)
  RealComplex.new(ExpS[x].Sum.Shift(n),ExpS[x].abs.Shift(n)) "}

これは matrix.rb を一寸修正すれば正方行列に対しても使えます。例えば 
Exp[[Matrix[0,-1],[1,0]]].printd(8) =>
Matrix[[0.54030230,-0.84147098],[0.84147098,0.54030230]]

Sqrt も複素数に対応しました。

Sqrt=Function.new{"|x| 
  a=(x<0)? Complex::I : 1
  s=Sequence.new([a],x){\"|n,x| a=self[n-1];(a+x/a)/2 \"}.Shift
  error=Sequence.new([],x,s){\"|n,x,s| y=s[n];(y-x/y).abs \"}
  RealComplex.new(s,error)
  "}

Sqrt[Complex::I].printd(16) =>
Complex(0.7071067811865475,0.7071067811865475)

以下は ε algorithm を含め今までの program の修正追加分です。

-----
class Function
  def *(y)
    x=self
    case y
    when Function
      Function.new({},x,y){"|z,x,y| x[y[z]]"}
    when Sequence
      Sequence.new([],x,y){"|n,x,y| x[y[n]] "}
    end 
  end
end

class Sequence
  def *(y)
    x=self
    case y
    when Numeric
      Sequence.new([],x,y){"|n,x,y| x[n]*y"}
    when Sequence
      Sequence.new([],x,y){"|n,x,y| x[n]*y[n]"}
    when Function
      Sequence.new([],x,y){"|n,x,y| x[y[n]]"}
    end
  end

  def alt
    AltSgn*self
  end

  def epsilon0
    x=self
    y=Sequence.new([],x){" |n,x| 1/(x[n+1]-x[n]) rescue nil "}
    Sequence.new([x,y]){" |n| a=self[n-2];b=self[n-1]
      Sequence.new([],a,b){\" |m,a,b|  a[m+1]+1/(b[m+1]-b[m]) \"} "}
  end

  def epsilon
    x=epsilon0
    Sequence.new([],x){" |n,x| x[2*n][0] "}
  end

  def Aitken
    x=self
    Sequence.new([],x){"|n,x| x[n+2]-((x[n+2]-x[n+1])**2)/(x[n+2]-2*x[n+1]+x[n
]) "}
  end

  def Aitken2
    Sequence.new([self]){"|n| self[n-1].Aitken "}
  end
end

AltSgn=Sequence.new([1,-1]){"|n| self[n-2]"}
Double=Function.new{"|x| 2*x"}
ExpS=Function.new{"|x| Power[x]/Fact "}
OddPower=Function.new{"|x| Power[x].Shift*Double "}

class Sequence2
  def initialize(s,e)
    case e
    when Sequence
      @seq=s
      @error=e
    when "D"
      @seq=s.Shift
      @error=@seq.Diff.abs
    when "E"
      @seq=s.epsilon.Shift
      @error=@seq.Diff.abs
    when "S"
      @seq=s.Sum.Shift
      @error=s.abs.Shift
    when "A"
      @seq=s.alt.Sum.Shift
      @error=s.abs.Shift
    when "W"
      @seq=Sequence.new([],s){"|n,s| s[n][0]"}
      @error=Sequence.new([],s){"|n,s| s[n][1]"}
    end
    @max=@seq+@error
    @min=@seq-@error
  end

  def add(a)
    x=self
    y=Sequence2(a)
    z=[x[0]+y[0],x.error[0]+y.error[0],0,0]
    s=Sequence.new([z],x,y){"|n,x,y| 
        i,j=self[n-1][2,2]
        if x.error[i]>y.error[j]
          i+=1
        else
          j+=1
        end
        [x[i]+y[j],x.error[i]+y.error[j],i,j]
        "}
    Sequence2.new(s,"W")
  end

  def +(y)
    case y
    when Sequence2
      add y
#     Sequence2.new(@seq+y.seq,@error+y.error)
    when Integer,Rational,Float
      Sequence2.new(@seq+y,@error)
    else
      x , y = a.coerce(self)
      x + y
    end
  end

  def -(y)
    case y
    when Sequence2
      add -y
#     Sequence2.new(@seq-y.seq,@error+y.error)
    when Integer,Rational,Float
      Sequence2.new(@seq-y,@error)
    else
      x , y = a.coerce(self)
      x - y
    end
  end

  def approx(eps)
    z=self[find_error{|e,i| e <= eps }[1]]
    return z if z.type!=Sequence2
    z.approx(eps)
  end

  def printd(n=16,r=10)
    approx(r**(-n)).printd(n,r)
  end
end
-----