\documentclass[12pt]{jarticle}
\topmargin -0.5in
\oddsidemargin -0in
\evensidemargin -0in
\textheight 9.5in
\textwidth 6in
\IfFileExists{my.sty}{\usepackage{my}}{}
\IfFileExists{graphicx.sty}{\usepackage{graphicx}}{}
\IfFileExists{epsfig.sty}{\usepackage{epsfig}}{}
\title{Risa/Asir の新しい分散表現多項式パッケージ}
\author{野呂 正行 \\ (神戸大理)}
\date{}
\begin{document}
\maketitle
\def\gr{Gr\"obner 基底}
\def\st{\, s.t. \,}
\def\noi{\noindent} 
\def\ve{\vfill\eject} 

\section{開発の経緯}
Risa/Asir においては、多項式は再帰表現または分散表現により保持される。
後者はグレブナー基底関連計算における基本的なデータ形式であり、
単項式を表す構造体 {\tt oMP} の linked list である。

\vskip 5mm
\begin{tabular}{cc}
\begin{minipage}{.5\hsize}
\begin{verbatim}
typedef struct oMP {
    struct oDL *dl;
    P c;
    struct oMP *next;
} *MP;
\end{verbatim}
\end{minipage}
&
\begin{minipage}{.5\hsize}
\begin{verbatim}
typedef struct oDL {
    int td;
    int d[1];
} *DL;
\end{verbatim}
\end{minipage}
\end{tabular}
\vskip 5mm

{\tt oMP} の係数は {\tt c} である。
{\tt oDL} のメンバー {\tt d}
は単項式の指数ベクトルを表しており、実際には変数の個数分の長さの
配列がセットされる。再帰表現された多項式は分散表現に変換され、
Buchberger, $F_4$, あるいは change of ordering などのアルゴリズム
ドライバにより処理される。

Risa/Asir のグレブナー基底計算においては、ペアの選択戦略、斉次化、モジュ
ラ計算、効率的容量除去などさまざまな効率化の工夫を採り入れることにより、
有理数体上での計算効率に関しては一定の評価を得てきたが、有限体上
での計算においては、Singular の最近の版と比較すると大きく性能が劣っていた。
また、有理数体上においても、多倍長演算に {\tt gmp} を使用している
Singular などのシステムでは、近年とみに高速化した
{\tt gmp} の性能と、Risa/Asir で
使用している自主開発の多倍長演算機能との性能差により、必ずしも Risa/Asir の
優位性が主張できなくなってきた。
一方で、PC に搭載できるメモリ量も数 GB に達し、CPU もどんどん高速化し、
グレブナー基底計算の応用範囲はどんどん大きくなっている。そこで、
これまでのさまざまな経験および、実装に関する最近の知見をもとに、
できる限り高速な分散表現多項式計算およびグレブナー基底計算を実現する
パッケージ {\bf nd} (New Distributed polynomial package) を新規に書くことにした。

\section{効率化の工夫}

Buchberger アルゴリズムに関しては、
Gebauer-Moeller の useless pair detection、sugar strategy などに
より、アルゴリズム的にはある程度固まったが、最近になっていくつか
実装に関する提案がなされた。今回の実装に採り入れたものについて
説明する. 

\begin{enumerate}
\item geobucket

これは、多項式の加算を効率化するための方法であり, 
\cite{Geo} で提案され, Macaulay2, Singular など多くのシステムで
採用され, 実際に効果があることが実証されている.
正規化計算では, 数多くの多項式の加算が行われるが, 非常に
項数の多い多項式に, 項数の比較的少ない多項式を繰り返し足すような場合,
項どうしの比較演算のコストが大変大きくなる. geobucket とは, 多項式を
要素とする配列 $g$ であって, 適当な整数 $b$ (例えば 2) に対し, $g[i]$
の多項式が高々 $b^i$ の項数を持つようなものである. $g$ に, 項数 $l$
($b^{i-1} < l \le b^i$) の多項式 $p$ を足す場合, まず $g[i]+p$ を
計算する. これの項数が $b^i$ 以下ならそのまま $g[i]$ を置き換え,
$b^i$ より大きければ $g[i+1]$ に加える, という操作を geobucket の
条件が満たされるまで続ける. これにより, 和に現われる多項式の項の総数を
$N$ とするとき, $O(N\log N)$ のコストで多項式の和が計算できる.


\item 可変長指数ベクトル

{\tt oDL} のメンバーでは, 単項式の変数の各指数を 32 bit 固定で表現して
いたが, 多くの場合これは過分であり, 結果として, 多項式のもつ情報量
よりも余分にメモリが必要となっていた. これに対して, 必要最小限の
bit 長を指数に割り当てておき, あふれが生じる際に, サイズを変更して
多項式を作りなおすというのがこの方法である. これは \cite{Singular}
で提案されている方法である.

\item 配列による多項式の保持

Buchberger アルゴリズムにおける基本操作は, $f-mg$ ($f$, $g$ は多項式,
$m$ は単項式) である. これをさらに $mg$ を作る操作と, 多項式の和を
作る操作に分解して考える. 後者は geobucket により高速化が可能である.
前者に関しては, どうしようもなさそうにも思えるが, この演算は, $g$ の
表現方法により計算効率が大きく異なる. 結論を言えば, 
$g$ が 単項式の linked list で表現されているより, 単項式 (これ自身
配列である) の配列として表現されているほうが, $mg$ の計算が高速
である. $g$ は, すでに計算された中間基底なので, 配列として保持すること
に問題はない. ただし, $mg$ は, 多項式加算にまわるので, linked list
として表示されていないと都合が悪い. 以上により, 多項式は状況に応じて
linked list と配列の 2 つの表現をとることになった.

\item 関数のインライン化, unrolling

開発が進むについて, ボトルネックとなる部分が次第に低レベルな部分に
なっていった. 特に問題となるのが, 単項式の指数ベクトルに対する操作
である. すなわち, 指数ベクトルの和, 差, 比較, divisibility などである.
これらの部分に関しては, インライン化, および unrolling の是非を
個別に実験により判断した. 

\item reducer のサーチのハッシュ化

項 $t$ を割り切る頭項をもつ中間基底 (reducer) $g_i$ のサーチも,
他部分の効率化が進むにつれ, そのコストが問題になってきた.
reducer としては, 経験上,
$i$ の小さい順から探して, $t$ を割り切る最初のものを用いるのが
よいとされる (例外もあるが). このため, $t$ の reducer はあれば一意
にきまる. 
$t$ の reducer $g_t$ が見つかったら, 
$t$ のハッシュ値 $h_t$ を計算して, ハッシュテーブルの $h_t$ の位置に,
$(t,g_t)$ を登録する. $t$ の reducer を探す際には, $h_t$ の位置
に登録されたデータから, $t$ の reducer を探して, もしあればそれを
用いればよい. 

\item 斉次の場合の効率化

一般には, 新たに生成された中間基底で, 既存の中間基底の正規化は行わないが,
入力が斉次の場合には, ある(weight つき)全次数の処理が終った時点で
その次数の中間基底どうしで inter reduction を行う. この場合, 頭項は
変化しないので, criteria への影響はなく, また, 低い全次数から順に
中間基底を生成していれば, 既に, 現次数までの簡約グレブナー基底の
すべての要素が得られているので, これまでに 0 に簡約された S-poly は
やはり新しい基底でも 0 に簡約される. この処理を行うことにより,
以降の計算が簡約基底により正規化されることになり, 正規化が効率化
されることが期待できる. 

\item メモリ管理

計算途中, さまざまな大きさの領域が繰り返し必要となる. 特に多く必要とさ
れるいくつかの構造体用領域は, garbage collector (GC) で得たものを自前
でフリーリスト管理している. これは, GC による allocation, collection 
が一定のコストを伴うためである.  この管理は nd パッケージ内で閉じてお
り, かつフリーリストの root を 0 にしておけば, いずれ GC により回収さ
れる.
\end{enumerate}

\section{基本データ構造}

\vskip 5mm
\begin{tabular}{cc}
\begin{minipage}{.5\hsize}
\begin{verbatim}
typedef struct oND {
    struct oNM *body;
    int nv,len,sugar;
} *ND;
\end{verbatim}
\end{minipage}
&
\begin{minipage}{.5\hsize}
\begin{verbatim}
typedef struct oNDV {
    struct oNMV *body;
    int nv,len,sugar;
} *NDV;
\end{verbatim}
\end{minipage}
\end{tabular}
\vskip 5mm

これらは, いずれも分散表現多項式を保持するための構造体である.
{\tt ND} は linked list 形式の, {\tt NDV} は配列形式の多項式
を表す.

\vskip 5mm
\begin{tabular}{cc}
\begin{minipage}{.5\hsize}
\begin{verbatim}
typedef struct oNM {
    struct oNM *next;
    union oNDC c;
    UINT dl[1];
} *NM;
\end{verbatim}
\end{minipage}
&
\begin{minipage}{.5\hsize}
\begin{verbatim}
typedef struct oNMV {
    union oNDC c;
    UINT dl[1];
} *NMV;
\end{verbatim}
\end{minipage}
\end{tabular}
\vskip 5mm

これらは, 単項式を表すための構造体である.  {\tt dl} は単項式の指数ベク
トルを表しており、実際には, 構造体作成時点での指数のbit 長と変数の
個数に応じた長さの配列の大きさ分の領域が確保される.
{\tt NM} は linked list 形式の, {\tt NMV} は配列形式の多項式における
単項式を表す. {\tt NDV} は, {\tt oNMV} すなわち構造体そのものの
配列へのポインタを持つ. 

\vskip 5mm
\begin{tabular}{cc}
\begin{minipage}{.5\hsize}
\begin{verbatim}
typedef union oNDC {
    int m;
    Q z;
    P p;
} *NDC;
\end{verbatim}
\end{minipage}
&
\begin{minipage}{.5\hsize}
\begin{verbatim}
typedef struct oRHist {
    struct oRHist *next;
    int index;
    int sugar;
    UINT dl[1];
} *RHist;
\end{verbatim}
\end{minipage}
\end{tabular}
\vskip 5mm

{\tt NDC} は係数を保持するための汎用の共用体である. 
{\tt m} は, 位数が 1 ワードで収まる有限体の元を保持するための
メンバーである. 
{\tt RHist} は reducer の履歴をハッシュテーブルに登録するための構造体である.
各エントリは, {\tt RHist} のリストとして登録される. 
\section{各部の詳細}

\subsection{ドライバ}

Buchberger アルゴリズムのドライバは, {\tt nd\_gb} と {\tt nd\_gb\_trace}
の二つがある. {\tt nd\_gb} は, 任意の係数体上で, sugar ストラテジーつきの
Buchberger アルゴリズムを実行するためのものである. 
ここでは, 

\begin{enumerate}
\item S-pair リストのメンテナンス
\item S-pair の取り出し, 正規化計算の呼び出し
\item 正規形の, content 除去, {\tt NDV} への変換
\item 指数にあふれが出た場合の, 中間基底の作りなおし
\end{enumerate}

などが行われる. {\tt nd\_gb\_trace} は, 有理数体, 有理関数体上のグレブナー
基底計算を trace アルゴリズムにより行うためのものであり, 上記の仕事に
加え, 結果をチェックする関数の呼び出し, homogenization, dehomogenization
も行われる.

さらに, 現状では有限体上のみであるが, $F_4$ ドライバ {\tt nd\_f4} も
実装した. S-pair, 中間基底の扱いに関しては {\tt nd\_gb} と同様である.
symbolic preprocessing は, 専用の geobucket が実装されている.
$F_4$ の核心である, 複数の S-pair から, reducer をまとめて
行列として掃き出す作業を行うまえに, 各 reducer により S-poly を
正規化している. この操作を行うために, 各 reducer を, 圧縮ベクトル
形式に変換しておき, 正規化される側の S-poly は非圧縮のベクトル
形式として正規化を行う. 最後に, 残った部分を集めて行列とし,
掃き出しを行っている. これらにより, できる限り使用メモリ量を押えて
いる.

\subsection{指数ベクトルの変更}

指数ベクトルの変更は, 指数の和であふれが生じたときに必要となる.
これが起こり得るのは, S-poly の計算と, 正規形の計算における,
単項式と多項式の積の計算においてである. この中での, 単項式
どうしの積の計算のたびにチェックするのは非効率的なので,
各中間基底に対し, 各変数に対する指数の最大値を記録しておき,
そのベクトルとの和があふれを起こす場合に作りなおしをしている.

\subsection{中間基底の demand loading}

{\tt dp} 系で提供されているのと同様に,
nd においても, 中間基底をディスク上の指定されたディレクトリに
置くことができる. 指定方法は {\tt dp} 系と同様 {\tt dp\_gr\_flags()}
で指定する. ファイルは {\tt dp} 系と同様の形式なので, {\tt bload()}
で読むことができる.

\subsection{有理数体上での content 除去}

正規化計算途中での content 除去は, 常に行われる. 現状では
頭係数が 2 倍 (固定) になったときに除去が行われる。

\section{性能}

一般に, 有限体上の計算の場合, {\tt nd\_gr} は {\tt dp\_gr\_mod\_main}
より数倍高速である. また, 問題にもよるが, {\tt nd\_f4} は
{\tt nd\_gr} の数倍程度高速な場合がある. おなじみの cyclic-$n$ で
比較すると表 \ref{tab:cyclic}のような結果を得る.

\begin{table}[hbtp] 
\begin{center}
\begin{tabular}{c||c|c|c|c}
 $n$        & {\tt nd\_gr} & Singular & {\tt nd\_f4} & {\tt dp\_gr\_mod\_main} \\ \hline
  7         &   5.1        &  5.0         & 1.8      & 17                      \\
  8         &   124        &  135         & 34       & 564                     \\
  9         &   27810      &  29725       & 3951     &    ---                  \\
\end{tabular}
\end{center}
\caption{$GF(31991)$ 上での DRL 順序グレブナー基底計算 (cyclic-$n$)}
\label{tab:cyclic}
\end{table}
このように, 少なくとも cyclic-$n$ では, nd の実装の効果が十分に現われている.
表 \ref{tab:janet} は, 種々のベンチマーク問題 \cite{janet} の計算時間を示す.
\begin{table}[hbtp] 
\begin{center}
\begin{tabular}{cc}
\begin{minipage}{.5\hsize}
\begin{tabular}{c||c|c|c}
       & {\tt nd\_gr} & Singular & {\tt nd\_f4} \\ \hline
dl & 5.9 & 4.9 &4.0 \\
eco10 & 7.1 & 10 &3.1 \\
eco11 & 63 & 106 &23 \\
eco12 & 507 & 1012 &198 \\
extcyc6 & 11 & 9.4 &4.1 \\
extcyc7 & 1813 & 1283 &447 \\
f855 & 3.6 & 3.4 &2.5 \\
filter9 & 0.28 & 0.80 &3.2 \\
hairer2 & 5.9 & 3.8 &4.5 \\
hairer3 & 11 & 35 &* \\
hcyclic7 & 6.5 & 4.8 &3.1 \\
hcyclic8 & 213 & 163 &82 \\
hf744 & 1.1 & 1.1 &1.6 \\
hf855 & 25 & 25 &17 \\
ilias13 & 11 & 8.4 &6.0\\
ilias\_k\_2 & 3.1 & 2.7 &1.1
\end{tabular}
\end{minipage}
&
\begin{minipage}{.5\hsize}
\begin{tabular}{c||c|c|c}
       & {\tt nd\_gr} & Singular & {\tt nd\_f4} \\ \hline
ilias\_k\_3 & 4.4 & 2.9 &1.2 \\
katsura10 & 285 & 218 &80 \\
katsura8 & 4.1 & 3.3 &1.3 \\
katsura9 & 35 & 29 &11 \\
noon7 & 4.4 & 1.8 &13 \\
noon8 & 35 & 18 &220 \\
pinchon1 & 3.6 & 1.0 &7.6 \\
rbpl & 1.0 & 0.89 &1.2 \\
redcyc7 & 3.5 & 3.3 &1.2 \\
redeco10 & 2.8 & 2.3 &1.3 \\
redeco11 & 24 & 18 &12 \\
redeco12 & 177 & 134 &74 \\
reimer6 & 11 & 32 &10 \\
reimer7 & 4000 & 4108 & 956 \\
virasoro & 1.8 & 1.4 & 0.65 
\end{tabular}
\end{minipage}
\end{tabular}

\end{center}
\caption{$GF(31991)$ 上での DRL 順序グレブナー基底計算}
\label{tab:janet}
\end{table}

有理数体上の計算の場合, 多項式や, 指数ベクトルの表現方法以外に, 途中あらわれる
係数の膨張の方が, 計算時間に大きく影響を与える場合が多い. この点では
{\tt nd\_gr\_trace} と {\tt dp\_gr\_main} とでは大差ないので割愛するが,
より悪くなることはない. 特に, weight を適切に設定することにより \cite{Kimura}, 
係数膨張に関してもより挙動のよい計算が可能となることに注意しておく.

\section{今後の予定}

{\tt dp} 系にあって nd にない機能として, 有理関数体係数のグレブナー基底
計算と, 有理数体上の $F_4$ 計算がある. なるべく早いうちにこれらを実装
したいと考えている. また, 
tangent cone アルゴリズムを用いた local ring での標準基底
計算も, reducer を探す関数を新たに用意することで対応可能と考えている.

\begin{thebibliography}{99}
\bibitem{Geo}
Yan, T., The Geobucket Data Structure for Polynomials.
Journal of Symbolic Computation, {\bf 25}, 3 (1998), 285-293.
\bibitem{Singular}
Sch\"onemann, H., Singular in a Framework for Polynomial Computations.
Joswig, M. and Takayama, N. (eds.), Algebra, Geometry, and Software Systems, 
Springer (2003), 163-176.
\bibitem{janet}
{\tt http://invo.jinr.ru/}. また {\tt http://www.symbolicdata.org}
にはさらに多くのベンチマーク問題がおいてある.
\bibitem{Kimura}
木村, 野呂, グレブナー基底計算のための weight 生成アルゴリズム.
本研究集会における発表 (2003).
\end{thebibliography}
\end{document}