Next: Assembly SIMD Instructions, Previous: Assembly Functional Units, Up: Assembly Coding [Index]
GMPにおける浮動小数点数の乗算は,性能の良くない整数乗算を使って実装されています。この際役立つのが64bitマシンにおける mpn_mul_1
, mpn_addmul_1
,mpn_submul_1
関数です。mpn_mul_basecase
関数は32bit, 64bitマシンの両方で活用されています。
IEEE 53bitの倍精度演算を使うと,整数乗算は高々53bitまでしか正しい結果が得られません。 64x64乗算を,8個の16x32->48 bitに分割すると簡便になります。6個の 21x32->53 bit乗算に分割することもできますが,この場合は最小の21bit分については符号ビットだけ使用することになります。
64bitマシンにおける mpn_mul_1
ファミリーは,処理開始時に一つのリムを3または4つに分割します。ループの内部では,多倍長オペランドを32bit単位で分割します。これらの符号なし32bit整数に分割されたものを浮動小数点数に効率よく変換できるかどうかはマシンによってかなり違いが出ます。読み込んだデータを整数に分割し,ゼロでパディングして64bit長に変換したり,あるいはメモリを経由して浮動小数点数に変換したりすることも可能になっています。
積の一部を64bitリムに変換するのが,符号付きの変換としては最良です。全ての値が2^53より小さいのであれば,符号付きでも符号なしでも同じ処理になりますが,多くのプロセッサでは符号なしの変換機能は持っていません。
以下に示す図は,64bitリムを用いてmpn_mul_1
やmpn_addmul_1
関数を使って16x32 bit積を行う時のものです。ひとまとまりのリムオペランドVは4つの16bit長に分割されます。こうしてつくられる複数リムオペランド Uはループの中で,二つの32bit部分を形成します。
+---+---+---+---+ |v48|v32|v16|v00| V operand +---+---+---+---+ +-------+---+---+ x | u32 | u00 | U operand (one limb) +---------------+ --------------------------------- +-----------+ | u00 x v00 | p00 48-bit products +-----------+ +-----------+ | u00 x v16 | p16 +-----------+ +-----------+ | u00 x v32 | p32 +-----------+ +-----------+ | u00 x v48 | p48 +-----------+ +-----------+ | u32 x v00 | r32 +-----------+ +-----------+ | u32 x v16 | r48 +-----------+ +-----------+ | u32 x v32 | r64 +-----------+ +-----------+ | u32 x v48 | r80 +-----------+
p32 と r32も,p48 and r48も浮動小数点加算可能です。 p00 と p16は,前の繰り返しで生成されたr64 と r80に加えられます。
それぞれのループで,4つの49bit長の結果を整数に変換し,下記のように並べられます。
|-----64bits----|-----64bits----| +------------+ | p00 + r64' | i00 +------------+ +------------+ | p16 + r80' | i16 +------------+ +------------+ | p32 + r32 | i32 +------------+ +------------+ | p48 + r48 | i48 +------------+
次に行うべきことは,これらを効率的に加算しつつ,桁上がりリムも加え,低位の64bitリムの結果と,高位の33bit桁上がりリムを生成していくかということです (i48は33bit分,高い桁数の方向に拡張されます)