多数桁計算用部品プログラム

　　　　　　2003/01/10　日立製作所 & 早稲田大学　後保範 ( Ushiro Yasunori )

1. 概要

　多数桁の加減算、乗除算等の多数桁の基本演算のFORTRANとCのプログラムを示す。
本プログラムは内容を理解してもらうことを目的とし、計算時間及び使用メモリ量の効率化
は目的としていない。本格的なプログラムは金田研究室から他のプログラムと合わせて
公開する予定。多数桁の乗算はFMT(高速剰余変換) を用いて作成している。
　多数桁の値はn個の倍精度実数で表わす固定小数点と、n+1個の倍精度実数で表わす
小数点の2形式がある。浮動小数点形式の最初の倍精度実数に指数部をもつ。
計算する値の基底は2進数でも10進数でも何でも良いが下記で最初に指定する。
　(1) FORTRANでの指定
　　　　COMMON /VBASE/BASE,NORID
　　　　BASE = 16.0**4　　! １要素が16進4桁の例 (10進5桁だと、10.0**5)
　　　　NORID = -1 　　! 各要素は符号付きでBASE/2以下にする。計算に利用。
　　　　NORID = 1　　　! 各要素を正符号でBASE以下にする。出力用の変換に利用。
　(2) Cでの指定
　　グローバル変数とし下記を指定する。10進5桁の例で示す。
　　　　double BASE=1.0e5;
　　　　int NORID=-1;
出力用の変換に利用するため、各要素を正符号にするには、下記を指定する。
　　　　NORID=1;
　乗算、除算及び平方根の計算では引数ERRに実数から整数への変換時の最大誤差を返す。
ERRが0.1を超える場合はBASEに指定する桁数を1つ小さくすると良い。ERRが0.5に近い値
の場合は結果が保証されない。通常はBASEは10進5桁か16進4桁が適当である。

2. 多数桁乗除算プログラム

　多数桁の乗算はFMT(高速剰余変換) を用いて作成している。演算結果は全て正規化している。

2.1 計算方法

　(1) 乗算
　　　FMT応用の「5. 複素FMTによる実数の畳込み演算」を利用する。
　　　本計算は(4)(a)のCFM4ルーチンで実現している。本計算の順変換をFMT、逆変換をFMT^-1
　　　とするとC=A･Bの多数桁乗算は下記のように計算される。
　　　　　Cの上位桁 = FMT^-1 ( FMT(A) × FMT(B) )の虚数部
　　　　　Cの上位桁 = FMT^-1 ( FMT(A) × FMT(B) )の実数部
　　　ここで、×は対応する要素単位の乗算を示す。
　(2) 除算
　　　x = 1/Aの計算を下記の反復計算で求める。
　　　　　x = x + x･( 1 - A･x )
　(3) 平方根
　　　x = SQRT(1/A)の計算を下記の反復計算で求める。
　　　　　x = x + x･( 1 - A･x² ) / 2

2.2 プログラム一覧

　(1) 乗算
　　(a) FMTCFM　：　浮動小数点の多数桁乗算 (C=A*B)
　　(b) MTCFM　　：　固定小数点の多数桁乗算 (C=A*B)
　(2) 除算
　　(a) FINV　　　：　浮動小数点の多数桁逆数計算 (C=1/A)
　　(b) FDIV　　　：　浮動小数点の多数桁除算 (C=A/B)
　(3) 平方根
　　(a) FCSQRT　：　浮動小数点の多数桁平方根 (入力は整数値)
　　(b) FSQRT　　：　浮動小数点の多数桁平方根 (C=SQRT(A) or 1/SQRT(A))
　(4) FMTルーチン
　　(a) CFMT4　　：　実数の畳込み演算用に変換した複素FFT
　　(b) CFMT　　　：　複素FMT
　　(c) CFMTTB　：　複素FMT用計算テーブル作成
　　(d) CFMTSB　：　複素FMTのステップ単位計算ルーチン

　ソースプログラムはFORTRAN とCを掲載する。

3. 多数桁加減算、正規化等プログラム

　演算結果は全て正規化している。

3.1 多数桁の正規化の方法

　(1) 固定小数点の正規化
　　　1要素(倍精度浮動小数点)にE(BASEの値、10⁵等)進数で与える。このとき、
　　　n要素の固定小数点の値fは下記のように表わされる。
　　　　　f = a₁･E^n-1 + …… + a_n-1･E + a_n
　　　計算中(NOIRD=-1)の正規化とは上記係数a_kを符号付きで下記の範囲に収める。
　　　　　| a_k | < E/2
　　　出力時の表示を考慮しNOIDR=1に設定して正規化すると係数a_kを下記にできる。
　　　　　0 < a_k < E
　　(2) 浮動小数点の正規化
　　　多数桁浮動小数点は先頭の要素(倍精度実数)が指数部で、その後にn要素の仮数部
　　　が続く構成である。指数部の値がαで仮数部がn要素の浮動小数点の値gは下記のよう
　　　に表わされる。
　　　　　g = E^α ･ ( a₁･E^-1 + …… + a_1-n･E + a_n･E^-n )
　　　浮動小数点の正規化は総ての係数a_kがゼロでない限り、先頭の係数 a₁がゼロ以外の
　　　数値になるよう仮数部をシフトし、指数部のαの値を調整する。

3.2 プログラム一覧

　(1) 加減算
　　(a) FADD　　：　浮動小数点の多数桁加減算 (C=A+B or A-B)
　　(b) ADD　　　：　固定小数点の多数桁加減算 (C=A+B or A-B)
　　(c) FCADD　：　浮動小数点の定数と多数桁の加減算 (C=k+B or k-B)
　(2) 正規化
　　(a) FNORM　：　浮動小数点の正規化（浮動小数点演算に利用している)
　　(b) NORM　　：　固定小数点の正規化（四則演算にに利用している)
　(3) 定数の乗除算
　　(a) FCMT　　：　多数桁と定数の乗算 (C=A*k)
　　(b) FCDIV　　：　固定小数点の正規化 (C=A/k)
　(4) その他
　　(a) COPY　　：　データのコピー
　　(b) SETV　　：　全要素に一定値を設定
　　(c) FVSET　：　定数を浮動小数点の多数桁の形式に変更
　　(d) LOG2F　：　整数に対して2のべき乗を求める

　ソースプログラムはFORTRAN とCを掲載する。

4. 時間計測プログラム

　ソースプログラムはFORTRAN を掲載する。
　Cはシステム提供関数(clock)を使用しているので不要。