2.8 その他(未分類)

2.8.1 tk_pfn.rkn

tk_pfn.rkn(F,X,Y,Xs,Ys,Ht,H)

:: Pfaffian 方程式に対する Runge-Kutta 法

return

リスト 独立変数と解の組

F, X, Y, Xs, Ys, Ht, H

F は Pfaffian 方程式の係数行列リスト. X は独立変数リスト. Y は従属変数リスト. Xs 独立変数の出発値リスト. Ys は出発時の従属変数の値リスト. Xt は停止する独立変数の値リスト. H は微少数.

• この関数は連立Pfaffian方程式 dY/d X[i] = F[i] Y を数値的に解く.
• 任意の holonomic system は Pfaffian 方程式に変換できる ([SST, Chap 1]).　変換には yang.rr パッケージを用いる.
• d F[i]/d X[j] + F[i] F[j] = d F[j]/d X[i] + F[j] F[i] = 0 が任意の i, j に対して成立していることが解が存在する必要十分条件である. この条件が成立しないときにこの関数を用いて解を計算してもその解は偽物である.
• X[i] が動く範囲は実数でないといけない.
• Xs[i] <= X[i] <= Xt[i] または Xt[i] <= X[i] <= Xs[i] である.
• 引数の与え方の例はソースコード (OpenXM/lib/asir-contrib/tk_pfn.rr )の tk_pfn.test1, tk_pfn.test2 を参照.
• 下の例の出力は X=(1,3) での値が Y=(-8,2,-6) であることを意味する.
• 参考. taka_runge_kutta.rr, yang.rr

[1355] import("tk_pfn.rr");
[1590]  tk_pfn.test1();
Value at (3,0.1)[8.99,6,-0.2]
Value at (1,3)[-8,2,-6]
[[[1,3],-8,2,-6],
 [[1,2.9],-7.41,2,-5.8],
 --- snip ---
 [[3,0.1],8.99,6,-0.2]]

ChangeLog

• この関数は 2009-12 から 2010-01 にかけて最初の版が書かれた.
• OpenXM/src/asir-contrib/packages/src/tk_pfn.rr 1.1, 1.2
• このモジュールの前身は tk_pf2.rr である. これは独立変数が 2 個の場合である.

2.8.2 tk_pfn.graph

tk_pfn.graph(Pf,Dom,Iv,Step)

:: 2変数 Pfaffian 方程式を Runge-Kutta 法で解いてグラフ表示する.

return

リスト リストの要素は以下の形式 [[xの値,yの値],Y_1の値,Y_2の値, ...]. [xの値, yの値] は [0,0],[0.2,0],[0.4,0], ... , [0,0.2],[0.2,0.2], ... のように y の値が外側ループ, x の値が内側ループの形式で増える.

Pf, Dom, Iv, Step

Pf は Pfaffian 方程式の係数行列リスト. 独立変数は x, y で固定. Dom リスト. 解くべき領域. Iv リスト. 領域の左端での初期値. Step 刻み幅.

• tk_pf2.rr, mt_graph.rr を import しておく必要がある.
• この関数は連立Pfaffian方程式 dY/dx = Pf[0] Y, dY/dy = Pf[1] Y を数値的に解いてグラフ表示する.
• Dom は [[xmin,xmax],[ymin,ymax]] の形式.
• 例はソースコード (OpenXM/lib/asir-contrib/tk_pfn.rr )の tk_pfn.testgraph1(), tk_pfn.testgraph2() を参照.
• option としては fit=1 がある. Z軸を適宜調整する.
• Dom の端はグラフ表示の時に一部カットされるので注意.

[1355] import("tk_pf2.rr"); import("mt_graph.rr"); import("tk_pfn.rr");
[1590] tk_pfn.testgraph1();

ここで testgraph1() は以下のとおり. 
def testgraph1() {
  /* tk_bess2.bess2pf(1/2); */
  Pf=  [[[ 0, (1)/(x), 0 ],
         [ -x, (2*x^2+1)/(x), -2*x ],
         [ -y, 0, 0 ]],
        [[ 0, 0, (1)/(y) ],
         [ -x, 0, 0 ],
         [ -x, (1/2)/(x), (-1/2)/(y) ]]];
  /* tk_bess2.bess2Iv(1/2,[0.5,1.5]); */
  Iv = [0.105994,-0.651603,-0.760628];
  Dom=[[0.5,1.5],[1.5,9]];
  Step = 0.5;
  return tk_pfn.graph(Pf,Dom,Iv,Step | fit=1);  
}

ChangeLog

• この関数は 2010-08 に最初の版が書かれた.
• OpenXM/src/asir-contrib/packages/src/tk_pfn.rr 1.8

2.8.3 tk_rk.runge_kutta_4

tk_rk.runge_kutta_4(Eq,X,Y,X0,Y0,Terminal,Step)

:: 4次の Runge-Kutta 法による微分方程式の数値近似解

return

リスト リストの要素は以下の形式 [Xの値,Y_1の値,Y_2の値, ...]. X の値は減っていく. よってリストの先頭が Terminal 付近での Y の値.

Eq, X, Y, Step

Eq は 方程式の右辺. リスト. Y[0]’=Eq[0], Y[1]’=Eq[1], ... である. X 独立変数名. Y リスト. 従属変数のリスト. Step 刻み幅.

X0, Y0, Terminal

X0 出発点の X の値. Y0 出発点での Y の初期値. Terminal X の終着点.

• taka_runge_kutta.rr を import しておく必要がある.
• この関数は連立常微分方程式 Y[0]’=Eq[0], Y[1]’=Eq[1], ... を数値的に解く.
• 例はソースコード (OpenXM/lib/asir-contrib/src/taka_runge_kutta.rr )の tk_rk.test4() を参照.

[1355] import("taka_runge_kutta.rr");
[1590] tk_rk.test4();

ここで test4() は以下のとおり. 振動の方程式, y0'=y1, y1'=-y0 (y0''+y0=0). 答は y0=cos(x)
taka_plot_auto は下方向で y が正.

def test4() {
  A=runge_kutta_4([y1,-y0],x,[y0,y1],0,[1,0],3.14*2,0.1);
  taka_plot_auto(A);
  return(A);
}

ChangeLog

• この関数は 2000 年代の前半に最初の版が書かれた. 2010年 Pfaffian の数値解析の為に再度整備
• OpenXM/src/asir-contrib/packages/src/taka_runge_kutta.rr 1.17

2.8.4 tk_rk.runge_kutta_4_linear

tk_rk.runge_kutta_4_linear(P,X,Y,X0,Y0,Terminal,Step)

:: 4次の Runge-Kutta 法による微分方程式の数値近似解. 線形方程式専用.

return

リスト リストの要素は以下の形式 [Xの値,Y_1の値,Y_2の値, ...]. X の値は減っていく. よってリストの先頭が Terminal 付近での Y の値.

P, X, Y, Step

P は 方程式の右辺. リスト. Y’=P Y である. 従属変数 Y は不要. X 独立変数名. Y リスト. 従属変数のリスト. 従属変数は自動生成される. 使われていないので空リストでよい. Step 刻み幅.

X0, Y0, Terminal

X0 出発点の X の値. Y0 出発点での Y の初期値. Terminal X の終着点.

• taka_runge_kutta.rr を import しておく必要がある.

[1355] import("taka_runge_kutta.rr");
[1590] A=tk_rk.runge_kutta_4_linear([[0,1],[-1,0]],x,[ ], 0, [1,0], 3.14*2, 0.1);
[1591] taka_plot_auto(A);

振動の方程式, y0'=y1, y1'=-y0 (y0''+y0=0). 答は y0=cos(x) を解いている.
taka_plot_auto は下方向で y が正.

ChangeLog

• 2010年 Pfaffian の数値解析の為に再整備.
• OpenXM/src/asir-contrib/packages/src/taka_runge_kutta.rr 1.17

2.8.5 fj_simplify.simplify

fj_simplify.simplify(arg1)

:: arg1 を簡単化する.

return

多項式, 有理式 または quote

arg1

多項式 または 有理式

• この関数は Joel S. Cohen, Computer Algebra and Symbolic Computation, http://web.cs.du.edu/~jscohen/MathematicalMethods/index.htm に記述されている Automatic simplification algorithm と B.F.Caviness, R.J.Fateman, Simplification of Radical Expressions (1976) に記述されている radcan アルゴリズムの実装である.
• 複素多値関数としては (x*y)^a = x^a*x^b は一般には成立しないので, 結果を複素関数に使うときは 注意が必要である. (x^a = exp(a*log(x)) なので)

import("fj_simplify.rr");
[1434] fj_simplify.simplify((x^(1/2))^3);
((x)^(3/2))
[1435] fj_simplify.simplify((2^(1/2))^2);
2
[1436] fj_simplify.simplify((2+2^(1/2))^3);
14*((2)^(1/2))+20
[1437] fj_simplify.simplify(exp(x)*exp(-x+y));
((@e)^(y))

参照

@ref{quote}

ChangeLog

• Todo: exp 以外の特殊関数についてのsimplification の機能.
• この関数は 2010.01 に M.Fujimoto により最初の版が書かれた. OpenXM/src/asir-contrib/packages/src/fj_simplify.rr

2.8.6 tk_jack.zonal

tk_jack.zonel(p,n)

:: 分割 p に対する n 変数の zonal 多項式を出力する.

return

多項式. x_1, x_2, ... が変数.

p

数のリスト. p=[p0,p1,...] の時 p0>=p1>=...>0.

n

自然数

• この関数は Koev-Edelman による Pieri 型公式を用いた Jack symmetric function の計算 アルゴリズムの実装である. 詳しくは Wikipedia 英語版 Jack symmetric function の項を参照.
• zonal(P,N) = jack(P,N,2) である.

load("tk_jack.rr");
[1434] tk_jack.zonal([3,2,1],3);
(112*x_3*x_2^2+112*x_3^2*x_2)*x_1^3+(112*x_3*x_2^3+168*x_3^2*x_2^2+112*x_3^3*x_2)*x_1^2+(112*x_3^2*x_2^3+112*x_3^3*x_2^2)*x_1
[1435] tk_jack.zonal([1,1],3);
(2*x_2+2*x_3)*x_1+2*x_3*x_2
[1436]  tk_jack.jack([1,1],3,2);
(2*x_2+2*x_3)*x_1+2*x_3*x_2

参照

ChangeLog

• この関数は wishart 分布に従う対称行列の第一固有値が x 以下である確率の計算を holonomic gradient method でやるためにその初期値を計算する C のプログラムが 必要であった. それを debug するためにとりあえず書いたもの.
• 最適化をまだまださぼってる.

2.8.7 ot_hgm_ahg.cbase

cbase(A)

:: A で定義される A-超幾何方程式系の Pfaffian の基底を求める.

return

Pfaffianの基底(微分作用素のモノミアル)のリスト

A

整数を成分とする行列 (maximal rank のもの)を表すリスト

• A-超幾何イデアルの Q(x)[dx] における標準基底は Pfaffian の基底となるが, 逆はかならずしも真ではない. 個数はもちろん同じである.
• アルゴリズムは T.Hibi, K.Nishiyama, N.Takayama, Pfaffian Systems of A-Hypergeometric Equations I, Bases of Twisted Cohomology Groups, arxiv:1212.6103 による. さらにパラメータ b を数に特殊化する確率算法を用いている.

[2190] import("ot_hgm_ahg.rr");
1
[2191] ot_hgm_ahg.cbase([[1,1,1,1],[0,1,2,3]]);
We use a probabilistic algorithm to determine the base.[dx2^2,dx3*dx2,dx3^2]

[dx3,dx4,1]

参照

@ref{get_mat2}

ChangeLog

• この関数は 2012 から 2014-春休みにかけてかかれた.
• version 1.1 以前の版は h-mle/A-hg/Prog (研究グループの共有フォルダ) にあり.

2.8.8 ot_hgm_ahg.get_mat2

get_mat2(A,W,Std,Mset)

:: A で定義される A-超幾何方程式系 H_A の Pfaffian の基底を求めるための Sylvester 法 を適用するための行列を生成する.

return

リスト

A

整数を成分とする行列 (maximal rank のもの)を表すリスト

W

リスト. toric ideal のグレブナー基底を計算するための weight vector. 多分なんでもいいはず.

Std

リスト. Pfaffian の基底. cbase(A) の出力を用いる.

Mset

Sylvester 型行列を作るための微分作用素のモノミアルのリスト.

• 出力を P に代入すると, P[0]*P[2] - P[1]*Std が modulo H_A で 0 となる. P[0] と P[2] を結合した行列が, sylvester 行列 (論文の記号での F’). P[2], Std が index モノミアルである. 論文での記号では P[2] は M_t, Std は S.
• アルゴリズムは K.Ohara, N.Takayama, Pfaffian Systems of A-Hypergeometric Equations II, Holonomic Gradient Method による. 論文の行列 F’.

[2190] import("ot_hgm_ahg.rr");
1
[2191] A=[[1,1,1,1],[0,1,2,3]]$
Std=ot_hgm_ahg.cbase(A)$
W=[[dx1,1,dx2,1,dx3,1,dx4,1]]$
Mset=[1,dx1,dx2,dx3,dx4]$
[2192] ot_hgm_ahg.get_mat2(A,W,Std,Mset);

省略

参照

@ref{cbase}

ChangeLog

• この関数は 2012 から 2014-春休みにかけてかかれた.
• version 1.1 以前の版は h-mle/A-hg/Prog (研究グループの共有フォルダ) にあり.
• ソース ot_hgm_ahg.rr の test3(), test3b(), test4(), test5(), test6(), test6c() 等に利用例がある.
• test3b() で Mset を一次式全部にしたものが, 論文の例.

2.8.9 ot_hgm_ahg.hgm_ahg_contiguity

hgm_ahg_contiguity(A,StdMon,Line,X0,InitVal,Start,End)

:: A で定義される A-超幾何方程式系のcontiguity relation を Sylvester matrix を用いて計算し, それを用いて超幾何関数の値を求める.

return

基底を超幾何関数に作用させたベクトルの値 F(End;X0) ??

A

整数を成分とする行列 (maximal rank のもの)を表すリスト.

StdMon

リスト. Pfaffian の基底を与える微分作用素のモノミアルのリスト.

Line

リスト [ContiDir,Beta,Z].

X0

リスト. x 変数の値.

InitVal

リスト. 基底を超幾何関数に作用させたベクトルの初期値 F(Start;X0)

Start

リスト. Z パラメータの初期値??

End

リスト. Z パラメータの終端値??

• Todo, この関数のインタフェースは変更される予定.
• X0 は有理数のリスト.
• ContiDir は End-Start と同じ方向.
• Beta. A超幾何関数の B パラメータの初期値 ??
• Z. ContiDir での一次元 contiguity を表現するための不定元の名前.
• ソース中の利用例. test_fd_conti(), test_c111_conti()
• アルゴリズムおよび利点は K.Ohara, N.Takayama, Pfaffian Systems of A-Hypergeometric Equations II, Holonomic Gradient Method 参照.

[2190] import("ot_hgm_ahg.rr");
1
[2191] ot_hgm_ahg.test_fd_conti();
(Todo, 引数がどうなるかの例を加える.)

参照

@ref{get_mat2} @ref{hgm_ahg_expected_value_contiguity} @ref{hgm_ahg}

ChangeLog

• この関数は 2014-07-11 に最初の版が 1.10版 ot_hgm_ahg.rr に commit された.
• インタフェースが更新された版は, 1.??版.

2.8.10 tk_hgpoly.optip

hgpoly.optip(A,B,W)

:: 整数計画問題をグレブナー基底を用いて解く.

return

リスト.

A

非負整数を成分とする行列 (maximal rank のもの)を表すリスト

B

非負整数を成分とするベクトルを表すリスト

W

非負整数を成分とするベクトルを表すリスト

• A U = B を満たす非負の整数ベクトル U の中で, 内積 W U を最小化する U を戻す.

[0] import("tk_hgpoly.rr");
[2191] tk_hgpoly.optip([[1,1,1,1],[0,1,2,3]],[20,40],[1,1,1,0]);
[6,1,0,13]

参照

@ref{feasible}

ChangeLog

• この関数は 2014-12-12 に commit された. 元版は h-mle/A-hg/Prog

2.8.11 tk_hgpoly.hgpoly

hgpoly.hgpoly(A,B)

:: A, B で定義される超幾何多項式を計算する.

return

リスト.

A

非負整数(todo, 再度確認)を成分とする行列 (maximal rank のもの)を表すリスト

B

非負整数を成分とするリスト.

• 戻り値を F とするとき, F[0] が超幾何多項式. 変数は x_1, x_2, ... F[1] は F[0] の分散表現多項式.

[0] import("tk_hgpoly.rr");
[2191] tk_hgpoly.hgpoly([[1,1,1,1],[0,1,2,3]],[2,2]);
[x_3*x_1+1/2*x_2^2,(1/2)*<<0,2,0,0>>+(1)*<<1,0,1,0>>]

ChangeLog

• この関数は 2014-12-12 に commit された.

2.8.12 tk_fd.abc2ahg

tk_fd.abc2ahg(A,B,C)

:: F_D(A, B,C) を解にもつ A-超幾何方程式系を求める.

return

リスト.

A

数

B

数のリスト

C

数

• 戻り値リストの第０成分は A-超幾何方程式系を定義する行列. 第１成分はA-超幾何方程式系のパラメータβ.

[2191] tk_fd.abc2ahg(-3,[-4,-5],3);

[[[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]],[11,5,4,5]]

参照

@ref{abc2marginal} @ref{marginal2abc}

ChangeLog

• この関数は 2014-12-13 に tk_fd.rr に追加された.

2.8.13 tk_fd.ahvec_abc

tk_fd.ahvec_abc(A,B,C,Y | all=1)

:: abc2marginal(A, B,C) を周辺和にもつ (2,m+1) 分割表全体についての正規化定数 Z, および Z の変数 Y[1][0], ..., Y[1][m] (2番目の行)についての偏微分を計算する.

return

リスト Ans

A

数

B

数のリスト. 長さは m.

C

数

Y

(2,m+1) 変数値をあらわすリストのリスト.

• A, B, C, に現れる数は整数を与える. Y の成分は有理数を与える.
• Ans[2]*Ans[1] が Z. Ans[2]*Ans[0][I] が Z の Y[1][I] についての偏微分.
• 論文 1.Y.Goto, Contiguity relations of Lauricella’s F_D revisited, arxiv:1412.3256 で導出されている contiguity relation を用いて計算する.

[2449] marginal2abc([3,12],[6,3,3,3]);
[-3,[-3,-3,-3],4]

[2450] tk_fd.ahvec_abc(-3,[-3,-3,-3],4,[[1,1/2,1/3,1/4],[1,1,1,1]]);
[[ 24041/1152 143551/11520 16973/1280 78827/5760 ],1/7776]

[2451] expectation_abc(-3,[-3,-3,-3],4,[[1,1/2,1/3,1/4],[1,1,1,1]]);
[721230/173593,430653/173593,458271/173593,67566/24799]

参照

@ref{expectation_abc}

ChangeLog

• この関数は 2014-夏に開発された.

2.8.14 pari, setbprec,setround,todouble,mpfr_gamma,mpfr_floor,mpfr_round,

pari(arg1,arg2,...)

:: MPFR で pari の関数を emulate するか ox_pari サーバーを呼び出す.

setbprec(arg1)

setround(arg1)

丸めの方法の指定. mpfr 準拠.

todouble(arg1)

:: bigfloat 型を double に変換する.

mpfr_gamma(arg1)

:: gamma 関数の計算.

mpfr_floor(arg1)

:: floor の計算.

mpfr_round(arg1)

:: 丸めの計算.

return

後述.

arg1, arg2

後述.

• 20150807 以降の asir では, bigfloat が pari ではなく mpfr を用いて計算される. Todo, pari は ox_pari を呼び出すが, まだ ox_pari に未実装の機能が多くある. 区間演算への対応は行っていない.
• pari(floor,arg1), pari(round,arg1), pari(gamma,arg1) は mpfr で pari を emulate しているので, 動作が異なる. 特に gamma は複素数の引数を受け取らない. pari の floor は桁数が足りなくなるとエラーで終了したが, この floor では setprec で指定した桁数以内なら正しく floor を戻す. Todo, ox_pari を指定して呼び出す方法.

[219] ctrl("bigfloat",1)$
[220] setprec(100)$
[221] pari(floor,1111111111111.1+1/10);
1111111111111

参照

@ref{pari}

ChangeLog

• これらの大変更は 2015-08-03 より 08-07 の asir 合宿でおこなわれた. まだ作業中. Todo, asirgui への対応. debug, ...
• 変更をうけたソースコードは asir2000/builtin/bfaux.c, asir2000/engine/bf.c, asir2000/builtin/parif.c 等多数.

2.8.15 xyz_pqr, syz_stu

xyz_pqr(arg1,arg2[,flag])

xyz_stu(arg1,arg2)

:: xyz に関する操作.

return

整数

arg1, arg2

整数

flag

0 または 1

• この項目は新しい関数の説明を書くためのテンプレートである. 消すな.
• xyz_pqr() は, arg1, arg2 を pqr する.
• flag が 0 でないとき, モジュラ計算を行う.
• xyz_stu() は stu アルゴリズムを用いる.

[219] xyz_pqr(1,2);
3
[220] xyz_pqr(1,2,1);
3
0
[221] xyz_stu(1,2);
3

参照

@ref{xyz_abc}

ChangeLog

• この関数は 2004-3-1 から 2004-3-14 にかけて アルゴリズム xyz (論文 http://www.afo.org/xyz.pdf ) を用いて書き直された. 変更をうけたソースコードは xxxyy.rr, ppp.c である.
• この関数は 2000 頃にはじめてのバージョンが書かれた. ソースは ppp.c である.

この文書は4月 24, 2024にtexi2html 5.0を用いて生成されました。