【内容一部公開】究極の精度を実現するために――近刊『有理算術演算—高精度数値計算のためのアルゴリズムとプログラミング—』
2023年9月下旬発行予定の新刊書籍、『有理算術演算—高精度数値計算のためのアルゴリズムとプログラミング—』のご紹介です。
同書の一部を、発行に先駆けて公開します。
＊＊＊
序文
ゴードン・ムーアは、1965年に集積回路あたりの部品数が毎年倍になると予測し、1968年にインテル社を創業した。1975年には次の10年を見据えて、2年ごとに2倍になるとの予測に修正した。予測はその後も維持され、「ムーアの法則」として広まり、計算機業界でもそれを目標に開発を進めた。この性能向上が維持されると、20年で千倍、40年では百万倍に達する。本書を執筆している時点の40年前は1982年で、数値計算の多くはメインフレームで実行されていた。当時の計算機の性能は、MIPS（Million Instructions Per Second）を基準に測られていた。当時を1MIPSとするなら、現在は百万MIPS、つまり1Tera Instructions Per Secondになる。しかし、この性能をそのまま実感できる人は少ないだろう。なぜなら、計算機の性能向上はクロック周期の短縮ではなく、マルチコア化に向かったため、十分に並列化されたプログラムを並列計算機で実行しないと、計算機のピーク性能に近づけないからだ。
数値計算は計算機科学の黎明期からの応用分野であり、数値シミュレーションの応用分野で多くのアプリケーションがFORTRANにより開発され、広く使用されるようになった。その中で、数値を扱う変数は、単精度と倍精度の浮動小数点数が主であった。この枠組みで数値シミュレーションは、計算機の能力向上の恩恵を十分に受けることができ、たとえば、数値線形代数では、非常に大きなサイズの行列を高速に処理できるようになった。40年前と比較すると、MIPS換算で百万倍に達する分野もあるかもしれない。
一方、高速化の対極にある高精度化の観点から見ると、この向上した計算能力はどのように生かされているだろうか。メーカーのFortranを中心に4倍精度がサポートされ、それを超える多倍精度計算も研究されている。この方向を追究していくと、究極の精度である有理算術演算（rational arithmetic）に到達するのだが、こちらはそれほど一般的な環境で試されているとは思えない。
有理算術演算は、原理的には小学校の算数で教えられる「分数計算」である。2つの自然数を、自動拡張する配列に格納し、それを分母と分子に使用し、さらに符号と組み合わせて有理数（分数）として保持する。四則演算は分数計算なので、演算のたびに分母と分子の桁数が増加するが、このとき自然数を格納する配列は自動拡張する。メモリと計算機パワーに対する要求は桁違いに増加するが、この計算方式なら、無理数が出てこない数式の計算では計算誤差は発生しない。計算量のオーダーが高いので、大きな行列に対して使用するには、計算機の能力は依然として不足する。現状では、小規模な問題を十進BASICの有理数モードで使用する教育者や、汎用プログラミング言語で倍精度浮動小数点演算と同居させて、数式処理システム上で有理算術演算を併用する研究者はいるだろうが、少数であろう。
BASICの文法では、数の表現は（PerlやPythonなどのスクリプト言語のように）1種類に限定されるので、有理数モードでの数（有理数）と、ほかの計算モードでの浮動小数点数を混在させたプログラムは書けない。筆者はC++言語で、有理数と浮動小数点数を同居させてプログラミング可能な「有理数計算プログラミング環境」を開発した。これは、通常のFortranなどのコンパイラ型のプログラミング言語が備えている数値（integer型やdouble型など）に加えて、有理数（rational型）、区間数（interval型）、可変精度の多倍精度浮動小数点数（mulprec型）、剰余数（モジュラー算術演算用のmodular型）などを備えて、これらの型に対する算術演算、論理演算、型変換などをサポートする。これによって、有理数解に対しては正確な解が求められ、無理数解に対してはこれを上下の有理数で挟み込む区間算術演算を用いて、必要な精度の解を挟む区間を求められる。プログラムは、計算内容が有理数の範囲か、無理数解を求めるものかを意識して書く必要がある。たとえば、係数行列と右辺ベクトルが有理数で与えられた連立1次方程式は、有理数の解ベクトルが誤差なしで得られる。対称行列が有理数で与えられた固有値問題では、特性多項式は有理数係数なので正確に求めて、特性方程式の解を区間演算で追究する。固有値が有理数の場合は正確な固有値を、無理数の場合は必要な精度まで区間を縮小反復するので、精度は保証しやすい。たとえば、浮動小数点演算でヒルベルト行列の固有値を求めると、（正の固有値が得られるはずなのに）負の固有値が得られることがある。このとき、行列が誤差を含んでいるためなのか、解法の問題なのかの特定は難しい。有理算術演算は部分的に適用することができるので、誤差を含んだ行列を有理数変換して、正確に解くことで、どちら側に問題があるかを特定できる。
「有理数計算プログラミング環境」によって、数値計算分野に、次に示すような新しい教育方法や問題分析の手法が開ける。
丸め誤差の理論を学ぶ前の学生を対象とした、数学教育
機械系の力学（弾性学や振動学）、電磁気学、量子力学などの数値シミュレーション分野の力学をカリキュラムにもたない数学科、数理科学科の学生を対象とした、HPC（High Performance Computing）プログラミング教育
数値シミュレーションで用いられる、浮動小数点演算による数値計算で発生した精度に関係する問題の分析
計算機援用証明（computer assisted proof、2つのアルゴリズムの計算結果の一致を確かめる）
ベンチマークの検収など、計算機システムの稼働確認
2番目の項目は、筆者自身の、教員の仕事で必要とされた。数理科学科の学生に「計算機がどうやって計算しているのか」を教えるためには、実習（とくに計算時間の測定）が有効である。学生はソフトウェア技術者として実務に就き、数十年間プログラミングをするかもしれない。その間、計算機の命令セットアーキテクチャの革新が1度や2度はあるだろう。このような問題意識に立つと、桁数の多い数値に対する四則演算プログラムを各自書いてみる演習は、「計算機がどうやって計算しているのか」を知るのによい入門になる。計算時間のかかる問題を解くことで、並列化の重要性も認識でき、コンパイル後にアセンブリ言語（機械語）が見えるコンパイラ型の言語で並列化することも実習できる。計算機システムに興味をもつことで、思考範囲が次世代計算機の命令セットアーキテクチャに及ぶ人材も育てたかった。40年後に、次の百万倍が達成されれば、有理算術演算は一般的な変数型として利用されるようになるだろう。
（後略）
＊＊＊
早稲田大学 寒川 光 (著)
筑波大学 高橋 大介(著)
【目次】
第I部 有理算術演算の構成
第1章 数の分類と「有理数計算プログラミング環境」の変数型
1.1 位取り記数法と基数変換
1.1.1 数の分類
1.1.2 基数変換
1.2 「有理数計算プログラミング環境」の有理数とその他の数
1.2.1 可変長非負整数longintと有理数rational
1.2.2 2進浮動小数点数と有理数
1.2.3 「有理数計算プログラミング環境」の変数型
1.2.4 「有理数計算プログラミング環境」での数値線形代数計算プログラム例
1.3 まとめ
第2章 有理算術演算の実装
2.1 プログラミング言語C++
2.1.1 C++小史
2.1.2 C++を採用したことによるエラー検出ルーチン
2.2 変数型と基本算術演算
2.2.1 可変長非負整数：longint数
2.2.2 有理数：rational
2.2.3 最大公約数（GCD）の計算
2.2.4 素因数分解の可能性
2.2.5 区間数：interval
2.2.6 多倍精度浮動小数点数：mulprec
2.2.7 53ビット精度の浮動小数点数：rough
2.2.8 mulprec型の数の平方根：SQR関数
2.2.9 プログラム例：桁数の超過
2.3 有理数の10進数表示
2.3.1 有理数rationalを表示する関数
2.3.2 区間数intervalを表示する関数
2.4 テンプレートメタプログラミング
2.4.1 Cのプリプロセッサによるマクロプログラミングの改良としてのテンプレート
2.4.2 再帰テンプレートによるループ構造のテンプレートプログラミング
2.4.3 CRTP概説
2.5 コンテナ（BLAS対応のベクトルと行列）
2.5.1 rtraits.hによる型の定義
2.5.2 rvector.hによるvectorクラス
2.5.3 rmatrix.hによるmatrixクラス
2.6 BLASアルゴリズム
2.6.1 rblasのメニュー
2.6.2 rblasの実装
2.6.3 prblasの実装
2.6.4 行列ベクトル積gemvの例
2.6.5 行列行列積gemmの例
2.6.6 並列化の例
2.7 まとめ
第3章 数値線形代数計算
3.1 連立1次方程式の直接解法
3.1.1 基本変換行列によるLU分解
3.1.2 LU分解プログラム
3.1.3 代入計算
3.1.4 連立1次方程式の解と固有値・固有ベクトルの関係
3.1.5 同次連立1次方程式の非自明解
3.1.6 LU分解と代入の性能
3.2 共役勾配法
3.2.1 共役勾配法
3.2.2 重複固有値をもつ問題
3.2.3 特異行列の場合
3.2.4 桁数削減による高速化
3.3 固有方程式の浮動小数点演算による解法
3.3.1 特性方程式
3.3.2 ヤコビ法
3.3.3 ギブンス法
3.3.4 ハウスホルダー法
3.3.5 対称3重対角行列の行列式の計算と2分法による固有値の計算
3.3.6 ギブンス法による多倍精度計算の精度と並列化の効果
3.3.7 ランチョス法
3.4 平方根の陰的定式化（3重対角化と線形最小2乗法）
3.4.1 共役勾配法を経由する3重対角行列の構成
3.4.2 線形最小2乗法
3.4.3 ギブンス法における平方根の陰的定式化
3.5 特性多項式と固有値の多重度
3.5.1 ヘッセンベルグ変換
3.5.2 フロベニウス変換
3.5.3 重根のある数値例
3.5.4 分割の細分化と固有値の多重度と重根に対応する固有ベクトル
3.5.5 特性多項式による3重対角行列の検定
3.6 特性多項式の分解
3.6.1 多項式の計算（ホーナー法）
3.6.2 多項式の格納と除算
3.6.3 根と係数の関係公式
3.6.4 因子探索プログラム
3.6.5 倍精度浮動小数点計算で作成された行列の固有値の多重度解析
3.6.6 近接根の分離
3.6.7 まとめ
3.7 スツルム関数列による固有値の計算
3.7.1 スツルム関数列の設定とスツルムカウントを数える関数
3.7.2 擬似乱数による行列の固有値
3.7.3 ヒルベルト行列の固有値
3.7.4 ヒルベルト行列の固有値計算の高速化
3.7.5 CT2D行列の近接根をもつ数値例
3.8 計算時間の振る舞いと有理算術演算の高速化
3.8.1 測定方法
3.8.2 ヘッセンベルグ変換の桁数と計算時間
3.9 まとめ
第4章 有理算術演算による無理数の計算
4.1 有理数・無理数と連分数
4.1.1 最良近似分数
4.1.2 一般の連分数展開とπの計算
4.1.3 自然数の平方根の連分数展開
4.2 平方根の計算
4.2.1 自然数（可変長非負整数）の開平演算
4.2.2 有理数の平方根（SQR関数とSQRint関数）
4.3 初等関数
4.3.1 ガウス-ルジャンドル法によるπの計算とgetpi関数
4.3.2 三角関数の引数還元
4.3.3 基本周期の三角関数
4.4 非線形方程式の求解
4.4.1 ratutilライブラリおよびmulutilライブラリの関数
4.4.2 2進数丸め関数roundratと挟み撃ち法bisectrf
4.4.3 フランク行列，ヒルベルト行列の固有値の精度検定と三角関数の精度検定
4.5 まとめ
第II部 有理算術演算の高速化
第5章 乗算の高速化
5.1 乗算アルゴリズム
5.1.1 ディジタル法
5.1.2 離散的フーリエ変換に基づくアルゴリズム
5.2 FFT乗算の実装
5.2.1 FFT乗算への分岐
5.2.2 lmulxとlmulfft関数
5.2.3 FFT関数
5.3 FFT乗算の計測
5.3.1 整数の行列による計測
5.3.2 GauLeg.cppによる計測
5.3.3 FHeig.cppによる計測
第6章 法算術演算による有理算術演算の高速化
6.1 法算術演算の基礎
6.1.1 中国剰余定理
6.1.2 拡張ユークリッド法と合同式の解法
6.1.3 連立合同式とその解法
6.2 法算術演算による有理算術演算の代替
6.2.1 有理算術演算を代替する法算術演算
6.2.2 有理数の復元
6.2.3 全剰余数の数え上げプログラム
6.2.4 法算術演算の作動確認プログラム
6.3 modularクラスによる法算術演算の実装と使用例
6.3.1 法の設定
6.3.2 不運な素数の数値例およびその対策
6.3.3 法算術演算によるLDL分解と不運な素数の検出と無効化
6.3.4 法算術演算によるCG法
6.4 法算術演算による特性多項式の計算の高速化
6.4.1 ヘッセンベルグ変換での中間式膨張と収縮のメカニズム
6.4.2 例題：charpoly
6.5 まとめ
第7章 まとめと将来展望
付録
付録A mempoolクラスによるlongint数の可変長配列による実装
A.1 可変長配列のサイズ
A.2 管理用配列
A.3 配列の定義と削除（コンストラクタとデストラクタ）
A.4 配列の拡張
A.5 マルチスレッド対応
A.6 コンストラクタ
A.6.1 コピーコンストラクタ
A.6.2 パラメータをとるコピーコンストラクタ
A.6.3 表示用コンストラクタ
付録B POSIXセマフォによるthread_managerの実装
付録C OpenMPによる並列化
C.1 3重対角化のtridiag.cpp（embarrassingly parallel）
C.2 スツルム関数列を使用したnumroots関数の並列化（配列の使用によるループ依存性の回避）
C.3 dvsrchプログラム（ループ変換）
C.3.1 プログラムの並列化
C.3.2 並列化プログラムの実行例
C.3.3 18分割データ
付録D ガウス超幾何関数による逆正接関数の連分数展開
付録E スツルム関数列
E.1 ユークリッド法で作るスツルム関数列
E.2 フーリエの定理によるスツルム関数列
E.3 スツルム関数列の係数の桁数
参考文献
索引
この記事が気に入ったらサポートをしてみませんか？