%comment $OpenXM: OpenXM/src/asir-contrib/packages/doc/n_wishartd/n_wishartd-en.texi,v 1.4 2021/03/29 03:12:40 takayama Exp $
%comment --- おまじない ---
\input texinfo
@iftex
@catcode`@#=6
@def@fref#1{@xrefX[#1,,@code{#1},,,]}
@def@b#1{{@bf@gt #1}}
@catcode`@#=@other
@end iftex
@overfullrule=0pt
@c -*-texinfo-*-
@comment %**start of header
@comment --- おまじない終り ---

@comment --- GNU info ファイルの名前 ---
@setfilename asir-contrib-n_wishartd

@comment --- タイトル ---
@settitle n_wishartd

@comment %**end of header
@comment %@setchapternewpage odd

@comment --- おまじない ---
@ifinfo
@macro fref{name}
@ref{\name\,,@code{\name\}}
@end macro
@end ifinfo

@iftex
@comment @finalout
@end iftex

@titlepage
@comment --- おまじない終り ---

@comment --- タイトル, バージョン, 著者名, 著作権表示 ---
@title n_wishartd
@subtitle n_wishartd User's Manual
@subtitle Edition 1.0
@subtitle Aug 2016

@author  by Masayuki Noro
@page
@vskip 0pt plus 1filll
Copyright @copyright{} Masayuki Noro
2016. All rights reserved.
@end titlepage

@comment --- おまじない ---
@synindex vr fn
@comment --- おまじない終り ---

@comment --- @node は GNU info, HTML 用 ---
@comment --- @node  の引数は node-name,  next,  previous,  up --- 
@node Top,, (dir), (dir)

@comment --- @menu は GNU info, HTML 用 ---
@comment --- chapter 名を正確に並べる ---
@menu
* n_wishartd.rr ::
* Index::
@end menu

@comment --- chapter の開始 ---
@comment --- 親 chapter 名を正確に ---
@node n_wishartd.rr ,,, Top
@chapter n_wishartd.rr 
@comment --- section 名を正確に並べる ---
@menu
* Restrction of matrix 1F1 on diagonal regions::
* Numerical comptation of restricted function::
* Differential operators with partial fraction coefficients::
* Experimental implementation of Runge-Kutta methods::
@end menu

This manual explains about @samp{n_wishartd.rr}, 
a package for computing a system of differential equations
satisfied by the matrix 1F1 on a diagonal region.
To use this package one has to load @samp{n_wishartd.rr}.
@example
[...] load("n_wishartd.rr");
@end example
@noindent
A prefix @code{n_wishartd.} is necessary to call the functions in this package. 

@comment --- section の開始 ---
@comment --- 書体指定について ---
@comment --- @code{} はタイプライタ体表示 ---
@comment --- @var{} は斜字体表示 ---
@comment --- @b{} はボールド表示 ---
@comment --- @samp{} はファイル名などの表示 ---

@node Restriction of matrix 1F1 on diagonal regions,,, n_wishartd.rr
@section Restriction of matrix 1F1 on diagonal regions 

@menu
* n_wishartd.diagpf::
* n_wishartd.message::
@end menu

@node n_wishartd.n_wishartd.diagpf,,, Restriction of matrix 1F1 on diagonal regions

@subsection @code{n_wishartd.diagpf}
@findex n_wishartd.diagpf

@table @t
@item n_wishartd.diagpf(@var{m},@var{blocks})
computes a system of PDEs satisfied by the @var{m} variate matrix 1F1 on a diagonal region specified
by @var{blocks}.
@end table

@table @var
@item return
A list @var{[E1,E2,...]} where each @var{Ei} is
a differential operator with partial fraction coefficients and it annihilates the restricted 1F1.

@item m
A natural number
@item vars
A list @var{[[s1,e1],[s2,e2],...]}.
@item options
See below.
@end table

@itemize @bullet
@item This function computes a system of PDEs satisfied by the @var{m} variate matrix 1F1 on a diagonal region specified by @var{blocks}.
@item Each component @var{[s,e]} in @var{blocks} denotes @var{ys=y(s+1)=...=ye}. The representative variable of this block is @var{ye}.
@item One has to specify @var{blocks} so that all the variables appear in it. In particular a block which contains only one variable
is specified by @var{[s,s]}.
@item It has not yet been proven that this function always succeeds. At least it is known that this function succeeds if the size of each block <= 36.
@item See @ref{Differential operators with partial fraction coefficients} for the format of the result.
@end itemize

@example
[2649] Z=n_wishartd.diagpf(5,[[1,3],[4,5]]); 
[
 [[[[-1,[]]],(1)*<<0,0,0,0,3,0>>],
  [[[-2,[[y1-y4,1]]],[-2,[[y4,1]]]],(1)*<<0,1,0,0,1,0>>],
  [[[9/2,[[y1-y4,1]]],[-3*c+11/2,[[y4,1]]],[3,[]]],(1)*<<0,0,0,0,2,0>>],
  ...
  [[[-6*a,[[y1-y4,1],[y4,1]]],[(4*c-10)*a,[[y4,2]]],[-4*a,[[y4,1]]]],
   (1)*<<0,0,0,0,0,0>>]],
 [[[[-1,[]]],(1)*<<0,4,0,0,0,0>>],

  [[[-6,[[y1-y4,1]]],[-6*c+10,[[y1,1]]],[6,[]]],(1)*<<0,3,0,0,0,0>>],
  [[[5,[[y1-y4,1]]],[-5,[[y1,1]]]],(1)*<<0,2,0,0,1,0>>],
  ...
  [[[21*a,[[y1-y4,2],[y1,1]]],[(36*c-87)*a,[[y1-y4,1],[y1,2]]],
   [-36*a,[[y1-y4,1],[y1,1]]],[(18*c^2-84*c+96)*a,[[y1,3]]],
   [-9*a^2+(-36*c+87)*a,[[y1,2]]],[18*a,[[y1,1]]]],(1)*<<0,0,0,0,0,0>>]]
]
@end example

@node n_wishartd.message,,, Restriction of matrix 1F1 on diagonal regions

@subsection @code{n_wishartd.message}
@findex n_wishartd.message

@table @t
@item n_wishartd.message(@var{onoff})
starts/stops displaying messages during computation.
@end table

@table @var
@item onoff
Start displaying messages if @var{onoff}=1.
Stop displaying messages if @var{onoff}=0.
@end table

@itemize @bullet
@item This function starts/stops displaying messages during computation.
The default value is set to 0.
@end itemize

@node Numerical comptation of restricted function,,, n_wishartd.rr
@section Numerical comptation of restricted function

@menu
* n_wishartd.prob_by_hgm::
* n_wishartd.prob_by_ps::
* n_wishartd.ps::
@end menu

@node n_wishartd.prob_by_hgm,,, Numerical comptation of restricted function
@subsection @code{n_wishartd.prob_by_hgm}
@findex n_wishartd.prob_by_hgm

@table @t
@item n_wishartd.prob_by_hgm(@var{m},@var{n},@var{[p1,p2,...]},@var{[s1,s2,...]},@var{t}[|@var{options}])
computes the value of the distribution function of the largest eigenvalue of a Wishart matrix.
@end table

@table @var
@item return
@item m
The number of variables.
@item n
The degrees of freedom.
@item [p1,p2,...,pk]
A list of the multiplicities of repeated eigenvalues.
@item [s1,s2,...,sk]
A list of repeated eigenvalues.
@end table

@itemize @bullet
@item 
Let @var{l1} be the largest eigenvalue of a Wishart matrix.
Let @var{Pr[l1<t]} be the distribution function of @var{l1}.
The function @code{n_wishartd.prob_by_hgm} computes the value of the distribution function by using HGM
for a covariance matrix which has repeated eigenvalues @var{si} with multiplicity @var{pi} (@var{i=1,...,k}).

@item This function repeats a Runge-Kutta method for the Pfaffian system by doubling the step size
until the relative error between the current result and the previous result is less than @var{eps},
The default value of @var{eps} is @var{10^(-4)}.
@item
If an option @var{eq} is not set, 
a system of PDES satisfied by 1F1 on the diagonal region specified by @var{[p1,p2,...]} is computed.
If an option @var{eq} is set, the list specified by @var{eq} is regarded as the correct system of PDEs.
@item If an option @var{eps} is set, the value is used as @var{eps}.
@item If an option @var{td} is set, the truncated power series solution for computing the initial vector
is computed up to the total degree specified by @var{td}. The default value is 100.
@item If an option @var{rk} is set, it is regarded as the order of a Runge-Kutta method. The default vaule is 5.
@item It is recommended to use this function only when @var{k<=2} where @var{k} is the number of diagonal blocks because of
the difficulty of the truncated power series solution and the difficulty of computation of the Pfaffian matrices.
@end itemize

@example
[...] n_wishartd.message(1)$
[...] P=n_wishartd.prob_by_hgm(10,100,[9,1],[1/100,1],100|eps=10^(-6));
...
[x0=,8/25]
Step=10000
[0]
[8.23700622458446e-17,8.23700622459772e-17]
...
Step=1280000
[0][100000][200000][300000]...[900000][1000000][1100000][1200000]
[0.516246820120598,0.516246820227214]
[log ratio=,4.84611265040128]

Step=2560000
[0][100000][200000][300000]...[2200000][2300000][2400000][2500000]
[0.516246912003845,0.516246912217004]
[log ratio=,4.93705929488356]
[diag,18.6292,pfaffian,1.09207,ps,41.0026,rk,213.929]
0.516246912217004
266.4sec + gc : 8.277sec(276.8sec)
@end example

@node n_wishartd.prob_by_ps,,, Numerical comptation of restricted function
@subsection @code{n_wishartd.prob_by_ps}
@findex n_wishartd.prob_by_ps

@table @t
@item n_wishartd.prrob_by_ps(@var{m},@var{n},@var{[p1,p2,...]},@var{[s1,s2,...]},@var{t}[|@var{options}])
computes the value of the distribution function of the largest eigenvalue of a Wishart matrix.
@end table

@table @var
@item return
@item m
The number of variables.
@item n
The degrees of freedom.
@item [p1,p2,...,pk]
A list of the multiplicities of repeated eigenvalues.
@item [s1,s2,...,sk]
A list of repeated eigenvalues.
@end table

@itemize @bullet
@item 
This function compute a truncated power series solution up to a total degree
where the relative error between the current value and the previous value at the desired point
is less than @var{eps}. The default value of @var{eps} is @var{10^(-4)}.
The value of the distribution function is computed by using this power series.
@item If an option @var{eps} is set, the value is used as @var{eps}.
@item
If an option @var{eq} is not set, 
a system of PDES satisfied by 1F1 on the diagonal region specified by @var{[p1,p2,...]} is computed.
If an option @var{eq} is set, the list specified by @var{eq} is regarded as the correct system of PDEs.
@item It is recommened to use this function when @var{t} is small.
@end itemize

@example
[...] Q=n_wishartd.prob_by_ps(10,100,[9,1],[1/100,1],1/2);
...
[I=,109,act,24.9016,actmul,0,gr,19.7852]
2.69026137621748e-165
61.69sec + gc : 2.06sec(64.23sec)
[...] R=n_wishartd.prob_by_hgm(10,100,[9,1],[1/100,1],1/2|td=50);
[diag,15.957,pfaffian,1.00006,ps,5.92437,rk,1.29208]
2.69026135182769e-165
23.07sec + gc : 1.136sec(24.25sec)
@end example

@node n_wishartd.ps,,, Numerical comptation of restricted function
@subsection @code{n_wishartd.ps}
@findex n_wishartd.ps

@table @t
@item n_wishartd.ps(@var{z},@var{v},@var{td})
computes a truncated power series solution up to the total degree @var{td}.
@end table

@table @var
@item return
A list of polynomial
@item z
A list of differential operators with partial fraction coefficients.
@item v
A list of variables.
@end table

@itemize @bullet
@item 
The result is a list @var{[p,pd]} where
@var{p} is a truncated power series solution up to the total degree @var{td}
and @var{pd} is the @var{td} homogeneous part of @var{p}.
@item @var{z} cannot contain parameters other than the variables in @var{v}.
@end itemize

@example
[...] Z=n_wishartd.diagpf(10,[[1,5],[6,10]])$
[...] Z0=subst(Z,a,(10+1)/2,c,(10+100+1)/2)$
[...] PS=n_wishartd.ps(Z0,[y1,y6],10)$
[...] PS[0];
197230789502743383953639/230438384724900975787223158176000*y1^10+
...
+(6738842542131976871672233/1011953706634779427957034268904320*y6^9
...+3932525/62890602*y6^2+1025/4181*y6+55/111)*y1
+197230789502743383953639/230438384724900975787223158176000*y6^10
+...+1395815/62890602*y6^3+3175/25086*y6^2+55/111*y6+1
@end example

@node Differential operators with partial fraction coefficients,,, n_wishartd.rr
@section Differential operators with partial fraction coefficients

@menu
* Representation of partial fractions::
* Representation of differential operators with partial fraction coefficients::
* Operations on differential operators with partial fraction coefficients::
@end menu

@node Representation of partial fractions,,, Differential operators with partial fraction coefficients
@subsection Representation of partial fractions

The coefficients of the PDE satisfied by the matrix 1F1 are written
as a sum of @var{1/yi} and @var{1/(yi-yj)} multiplied by constants.
Furthermore the result of diagonalization by l'Hopital rule
can also be written as a sum of partial fractions.

@itemize @bullet
@item 
A product @var{yi0^n0(yi1-yj1)^n1(yi2-yj2)^n2...(yik-yjk)^nk}
in the denominator of a fraction is represented as a list @var{[[yi0,n0],[yi1-yj1,n1],...,[yik-yjk,nk]]},
Where each @var{yi-yj} satisfies @var{i>j} and the factors are sorted according to an ordering.
@item 
Let @var{f} be a power sum as above and @var{c} a constant.
Then a monomial @var{c/f} is represented by a list は @var{[c,f]}.
@var{f=[]} means that the denominator is 1.
@item 
Finally @var{c1/f1+...+ck/fk} is represented as a list @var{[[c1,f1],...,[ck,fk]]},
where terms are sorted according to an ordering.
@item
We note that it is possible that a partial fraction is reduced to 0.
@end itemize

@node Representation of differential operators with partial fraction coefficients,,, Differential operators with partial fraction coefficients
@subsection Representation of differential operators with partial fraction coefficients

By using partial fractions explained in the previous section,
differential operators with partial fraction coefficients are represented.
Let @var{f1,...,fk} be partial fractions and @var{d1,...,dk} distributed monomials such that
@var{d1>...>dk}) with respected to the current monomial ordering.
Then a differential operator @var{f1*d1+...+fk*dk} is represented as a list @var{[f1,d1],...[fk,dk]]}.

@node Operations on differential operators with partial fraction coefficients,,, Differential operators with partial fraction coefficients
@subsection Operations on differential operators with partial fraction coefficients

@menu
* n_wishartd.wsetup::
* n_wishartd.addpf::
* n_wishartd.mulcpf::
* n_wishartd.mulpf::
* n_wishartd.muldpf::
@end menu

@node n_wishartd.wsetup,,, Operations on differential operators with partial fraction coefficients
@subsubsection @code{n_wishartd.wsetup}
@findex n_wishartd.wsetup

@table @t
@item n_wishartd.wsetup(@var{m})
@end table

@table @var
@item m
A natural number.
@end table

@itemize @bullet
@item This function sets a @var{m}-variate computational enviroment. The variables are @var{y0,y1,...,ym} and @var{dy0,...,dym},
where @var{y0, dy0} are dummy variables for intermediate computation.
@end itemize

@node n_wishartd.addpf,,, Operations on differential operators with partial fraction coefficients
@subsubsection @code{n_wishartd.addpf}
@findex n_wishartd.addpf
@table @t
@item n_wishartd.addpf(@var{p1},@var{p2})
@end table

@table @var
@item return
A differential operator with partial fraction coefficients.
@item p1, p2
Differential operators with partial fraction coefficients.
@end table

@itemize @bullet
@item This function computes the sum of differential operators @var{p1} and @var{p2}.
@end itemize

@node n_wishartd.mulcpf,,, Operations on differential operators with partial fraction coefficients
@subsubsection @code{n_wishartd.mulcpf}
@findex n_wishartd.mulcpf
@table @t
@item n_wishartd.mulcpf(@var{c},@var{p})
@end table

@table @var
@item return
A differential operator with partial fraction coefficients.
@item c
A partial fraction.
@item p
Differential operators with partial fraction coefficients.
@end table

@itemize @bullet
@item This function computes the product of a partial fraction @var{c} and a differential operator @var{p}.
@end itemize

@node n_wishartd.mulpf,,, Operations on differential operators with partial fraction coefficients
@subsubsection @code{n_wishartd.mulpf}
@findex n_wishartd.mulpf
@table @t
@item n_wishartd.mulpf(@var{p1},@var{p2})
@end table

@table @var
@item return
A differential operator with partial fraction coefficients.
@item p1, p2
Differential operators with partial fraction coefficients.
@end table

@itemize @bullet
@item This function computes the product of differential operators @var{p1} and @var{p2}.
@end itemize

@node n_wishartd.muldpf,,, Operations on differential operators with partial fraction coefficients
@subsubsection @code{n_wishartd.muldpf}
@findex n_wishartd.muldpf
@table @t
@item n_wishartd.muldpf(@var{y},@var{p})
@end table

@table @var
@item return
A differential operator with partial fraction coefficients.
@item y
A variable.
@item p
A differential operator with partial fraction coefficients.
@end table

@itemize @bullet
@item 
This function computes the product of the differential operator @var{dy} corresponding to a variable @var{y} and @var{p}.
@end itemize

@example
[...] n_wishartd.wsetup(4)$
[...] P=n_wishartd.wishartpf(4,1);
[[[[1,[]]],(1)*<<0,2,0,0,0>>],[[[1/2,[[y1-y2,1]]],[1/2,[[y1-y3,1]]],
...,[[[-a,[[y1,1]]]],(1)*<<0,0,0,0,0>>]]
[...] Q=n_wishartd.muldpf(y1,P);
[[[[1,[]]],(1)*<<0,3,0,0,0>>],[[[1/2,[[y1-y2,1]]],[1/2,[[y1-y3,1]]],
...,[[[a,[[y1,2]]]],(1)*<<0,0,0,0,0>>]]
@end example

@node Experimental implementation of Runge-Kutta methods ,,, n_wishartd.rr
@section Experimental implementation of Runge-Kutta methods 

@menu
* rk_ratmat::
@end menu

@node rk_ratmat,,, Experimental implementation of Runge-Kutta methods 

In the function @code{n_wishartd.ps_by_hgm}, after computing the Pfaffian matrices for
the sytem of PDEs on a diagonal region, it executes a built-in function
@code{rk_ratmat} which computes an approximate solution of the Pfaffian system
by Runge-Kutta method for a spcified step size.
This function is repeated until the result gets stabilized, by doubling the step size.
@code{rk_ratmat} can be used as a general-purpose Runge-Kutta driver and we explain how to use it.

@subsection @code{rk_ratmat}
@findex rk_ratmat

@table @t
@item rk_ratmat(@var{rk45},@var{num},@var{den},@var{x0},@var{x1},@var{s},@var{f0})
solves a system of linear ODEs with rational function coefficients.
@end table

@table @var
@item return
A list of real numbers.
@item rk45
4 or 5.
@item num
An array of constant matrices.
@item den
A polynomial.
@item x0, x1
Real numbers.
@item s
A natural number.
@item f0
A real vector.
@end table

@itemize @bullet
@item 
Let @var{k} be the size of an array @var{num}.
The function @code{rk_ratmat} solves  an initial value problem
@var{dF/dx = P(x)F}, @var{F(x0)=f0} for @var{P(x)=1/den(num[0]+num[1]x+...+num[k-1]x^(k-1))} by a Runge-Kutta method.
@item
@var{rk45} specifies the order of a Runge-Kutta method. Adaptive methods are not implemented.
@item
The step size is specified by @var{s}. The step width is @var{(x1-x0)/s}.
@item
If the size of @var{f0} is @var{n},  each component of @var{num} is a square matrix of size @var{n}.
@item
The result is a list of real numbers @var{[r1,...,rs]} of length @var{s}.
@var{ri} is the 0-th component of the solution vector after the step @var{i}.
Before going to the next step the solution vector is divided by @var{ri}.
Therefore the 0-th component of the final solution vector [var{F(x1)} is equal to @var{rs*r(s-1)*...*r1}.
@item Since the ODE is linear, each step of Runge-Kutta method is also linear.
This enables us to apply a normalization such that the 0-th
component of each intermediate solution vector is set to 1.  By
applying this normalization we expect that all the components of
intermediate solution vectors can be represented by the format of
double precision floating point number. 
If there exist some components which cannot be represented by the format
of double precision floating point number in the initial vector @var{f0}, we apply this normalization
to @var{f0}. After applying @code{rk_ratmat} we multiply the result for the normalized @var{f0} and the 0-th component
of the original @var{f0} to get the desired result.
@end itemize

@example
[...] F=ltov([sin(1/x),cos(1/x),sin(1/x^2),cos(1/x^2)]);
[ sin((1)/(x)) cos((1)/(x)) sin((1)/(x^2)) cos((1)/(x^2)) ]
[...] F0=map(eval,map(subst,F,x,1/10));
[ -0.54402111088937 -0.839071529076452 -0.506365641109759 0.862318872287684 ]
[...] N0=matrix(4,4,[[0,0,0,0],[0,0,0,0],[0,0,0,-2],[0,0,2,0]])$
[...] N1=matrix(4,4,[[0,-1,0,0],[1,0,0,0],[0,0,0,0],[0,0,0,0]])$
[...] N=ltov([N0,N1])$
[...] D=x^3$
[...] R=rk_ratmat(5,N,D,1/10,10,10^4,F0)$
[...] for(T=R,A=1;T!=[];T=cdr(T))A *=car(T)[1];
[...] A;
0.0998334
[...] F1=map(eval,map(subst,F,x,10));
[ 0.0998334166468282 0.995004165278026 0.00999983333416666 0.999950000416665 ]
@end example


@comment --- おまじない ---
@node Index,,, Top
@unnumbered Index
@printindex fn
@printindex cp
@iftex
@vfill @eject
@end iftex
@summarycontents
@contents
@bye
@comment --- おまじない終り ---