任意精度演算

任意精度演算（にんいせいどえんざん）^[1]とは、数値の精度を必要ならいくらでも伸ばしたりできるような演算システム（実際上は利用可能なメモリ容量に制限されるが）による演算である。

概要

多倍長整数（たばいちょうせいすう）などを内部処理に利用し、必要な桁数の浮動小数点計算を行う。固定長の整数や一般的な固定精度の浮動小数点方式は、ハードウェアで高速に処理できるのに対し、任意精度演算はソフトウェアで実装され、重い処理を必要とする。十進の0.1を2進で表現しようとする場合のように、有限の桁数では表現し切れない場合もあることから、2進でなく十進で処理するものや、有理数演算を併用したりもする。

多倍長演算（たばいちょうえんざん）^[2]とも言うが、プログラミング言語によっては、多倍長整数 (特に区別する場合は bigint などと言う) の名前が bignum であることもある。

最近のプログラミング言語の中には、多倍長整数を言語仕様でサポートしているものもあり、他の言語でも多倍長整数や任意精度の浮動小数点数を扱うライブラリが存在する。任意長の配列に格納するような実装になっている。

任意精度演算は、演算速度が重要視されない用途や大きな数についての正確な演算結果を必要とする場合に使われる。

数式処理システムとの違い

数式処理システムの記号計算（代数処理）では、たとえば a + b という式はそういう記号による式のまま表現し、たとえばその3倍は代数法則に従い 3a + 3b のように処理するので、任意の計算可能数を無限の精度で「表現」できるが、任意精度演算とは異なったものである。しかし、個々の数式処理システムの設計者の考え次第であるが、数式からシームレスに、あるいは明確なアクションを経て、任意精度演算や有理数演算に繋げられるものもある。

用途

典型的用途として数百から数千桁の整数の演算を必要とする、公開鍵暗号がある。また、人間中心のアプリケーションで桁数制限やオーバフローが不適切な場合にも使用される。また、固定精度演算の結果をチェックするのにも使え、方程式の係数（例えばガウス求積に出てくる $\pm {\sqrt {1/3}}$ など）の最善値を決定するのにも使える。

任意精度演算は、円周率などの基礎的数学定数を数百万桁以上計算して求め、その数字列の特性を解析するのに使ったり^[3]、解析的手法では解明するのが難しい関数（例えば、リーマンゼータ関数）の正確な振る舞いを調べるのに使ったりする。また、マンデルブロ集合などのフラクタル画像を極めて高い倍率で描く場合にも任意精度演算が必要となる。

任意精度演算は、固定精度演算では避けられない算術オーバーフローを防ぐために使うこともある。4桁の走行距離計が9999の次は0000に戻るように、固定精度の整数型では、数値が固定精度で表せる範囲を超えると循環するように最も小さい値に戻ってしまう。「循環」させずに「飽和」させる方式のプロセッサもあり、演算結果が表現できない数値になった場合に、表現可能な最も近い値に置換してしまう（例えば、16ビット符号なし整数なら、65535 に 1 を足しても 65535 のままとなる）。プロセッサによっては例外を発生させることもできる。必要なら例外をひきとって、ソフトウェアが任意精度演算に切り替えるということも可能である。

コンピュータの多くは整数として32ビットや64ビットを使うのが通例で、アプリケーションはオーバーフローを起こさないよう注意してプログラミングしなければならない。しかし、たとえば配列のサイズを扱っているのであれば、そのようなバグを踏むのはバイトの巨大な配列の時だけであり、しばしば忘れられている。例えば、二分探索のコードとして大抵の参考書でそのように書かれており、実際に世界中で広く使われていたコードで (L + R) / 2 という式で計算を行っていた。固定長の場合この計算の結果は L + R がオーバーフローした時に正しくないものとなる。しかし、もし任意長整数計算であったなら、この式のままで問題ないわけである。

LISP、REXX、Python、Perl、Ruby といったプログラミング言語では、整数演算で値の範囲が固定長の範囲を越えるものを、自動的に多倍長整数にフォールバックする（オプションの場合もある）。性能的には不利だが、演算結果がオーバーフローで不正（または例外）になる可能性を排除できる。また、ワードサイズの異なる機種間でも演算結果が常に同じになることを保証できる。プログラミング言語で多倍長整数をサポートすると、言語を単純化でき、様々なデータ型を用意する必要がないという利点もある。理論的には、最適化において数学的な式変形が可能になる、という利点もあるが、実際としてはフォールバックによる性能低下のほうが利いてくる。

実装上の問題

任意精度演算は、プロセッサのレジスタの大きさに合わせた数値型の演算より大幅に性能が劣る。固定精度演算はハードウェアでの実装が比較的容易であるのに対して、任意精度演算はメモリ管理などソフトウェアで実装しなければならない部分がある（多倍長ないし任意精度の演算のハードウェアサポートは過去広くみられる）。整数の除算や浮動小数点演算はサポートしていないプロセッサもあるが、ソフトウェアでそれらをサポートするとしても通常は1ワードまたは2ワードの数値型を使い、任意長の処理は行わない。

一般に除算では循環小数が発生することがある。任意精度演算ではどこかで打ち切るか、循環をなんらかの形で表現するか、有理数演算を併用する。

多倍長整数であれば任意長の整数演算を、任意精度浮動小数点であれば任意精度の浮動小数点の処理をおこなう。任意精度数の大きさは使用可能な記憶装置の容量で制限されるだけでなく、数値を表す配列のインデックスのサイズでも制限される。一般にインデックスは固定長である。

任意精度の数値の演算を効率的に行うためのアルゴリズムがいくつも考案されてきた。特に、桁数を $N$ としたとき、 $N$ が大きい場合の計算量を最小にするようなアルゴリズムが必要とされる。

加算や減算の最も単純なアルゴリズムは、単純に桁ごとに順次加算・減算していき、繰り上げや繰り下げを必要に応じて行うものである。この計算量は $O(N)$ である（ランダウの記号参照）。

乗算の場合、最も単純なアルゴリズムは筆算の計算方法をそのまま採用する手法で $O(N^{2})$ の計算量である。計算量が $O(N\log(N)\log(\log(N)))$ の乗算アルゴリズムとして、高速フーリエ変換に基づくショーンハーゲ・ストラッセン法（Schönhage-Strassen-Algorithmus）がある。また、カラツバ法では $O(n^{\log _{2}3})$ の計算量である。 $N$ があまり大きくない場合、カラツバ法の方が性能がよい（前者の性能が勝つのは少なくとも1万桁以上の場合）。

例

階乗の計算では、非常に素早く非常に大きな数を生成する。方程式などでは他の項との組合せで出現するので、評価の順序を工夫することで、全体の計算では多くの場合（テイラー級数など）それが問題になることはない。階乗の近似値が必要なら、スターリングの近似を使えばよい。正確な値が必要なら、整数型の限界をすぐに超えることが問題となる。浮動小数点数による近似であってもその最大値を超えるのは容易で、対策として階乗の対数による計算に置き換える方法が出てくる。

大きな階乗の正確な値を求めたい場合、特別なソフトウェアが必要となる。以下の擬似コードは、1、1*2、(1*2)*3、((1*2)*3)*4 と順に計算していく手法を使ったものである。

program Factorial ;
  type natural = range 1..integer.last of integer ;
  type radical = range 2..integer.last of integer ;
  type COLUMNS_TYPE = record
    const MAX_LENGTH : natural = 1000 ; (* 十分な桁数を上限とする *)
    values : array [1:MAX_LENGTH] of natural ;
    length : natural ;
    radix : radical ;
    end ;
  
  procedure factorial (const n : natural, var columns : COLUMNS_TYPE) ;
    begin
    columns.length := 1 ;
    columns.values [columns.length] := 1 ; (* 最初は乗法の単位元 1 を設定する *)
    var carry : integer = 0 ;              (* 下の位からの桁上り *)
    for var ci := 1 to columns.length do   (* 桁ごとにループ *)
      begin
      const c = columns.values [ci] * n + carry ;  (* この桁の数との乗算と下の位からの桁上りの和 *)
      columns.values [ci] := c mod columns.radix ; (* 積の下位桁 *)
      carry := c div columns.radix ;               (* 積の上位桁、上の位への桁上り *)
      end ;
    while 0 < carry do
      begin
      if COLUMNS_TYPE.MAX_LENGTH <= columns.length then
        raise OverflowError ;
      columns.length := columns.length + 1 ;                       (* 桁を増やす *)
      columns.values [columns.length] := carry mod columns.radix ; (* 実際に格納 *)
      carry := carry div columns.radix ;                           (* 次の位への桁上り *)
      (* n が基数より大きければ、もう1桁必要になることがある。 *)
      end
    end ;
  
  procedure writeColumns (var column : COLUMNS_TYPE) ;
    begin
    if 36 < columns.radix then raise OutOfRangeError ;
    const DIGITS : array of character = ['0'..'9','A'..'Z'] ; (* 36進法まで対応可能 *)
    for var ci := columns.length downto 1 do
      write (DIGITS [columns.values [ci]]) ; (* 大きい桁から表示 *)
    end ;
  
  begin
  for var n := 1 to 35 do
    begin
    var columns : COLUMNS_TYPE ;
    columns.radix := 10 ; (* 10進法 *)
    factorial (n, columns) ;
    writeColumns (columns) ; writeln (' = !', n)
    end
  end.

この例で、詳細を見てみよう。最も重要なのは、大きな数の表現方法の選択である。この場合、階乗の計算なので整数だけでよく、固定小数点方式で事足りる。基数の冪は0から始まって大きくなっていくので、配列の先頭が0で、後ろの方ほど基数の冪が大きくなるという表現が便利である。配列のインデックスの始点や大きさの指定方法は言語によって異なるが（0から始まる言語と1から始まる言語がある）、一般に計算の必要条件には関係しないので、ここでは単純化するために配列の始点を0ではなく1としている。数字配列のインデックスがその桁の基数の冪と対応するという性質は、ここでは利用していない。

次に重要な決定は、基数を10としたことである。この場合、様々なことを考慮しなければならない。一時変数 c は、1つの桁の乗算結果と前の桁の積からの繰り上がりを加えたものを保持しなければならない。基数が10の場合、16ビット整数であれば32767まで表現できるので十分である。ただし、この例では n 自身が10を基数とする配列になっていないという点で若干ごまかしている。結果としてこの例は n ＞ 3200 などと設定するとうまく動作できなくなる。これがこの擬似コードの暗黙の限界である。一般には n も大きな数として配列で表す必要がある。また、このようなごまかしをしたせいで、一桁の乗算であっても n 全体をかけているため、繰り上がりは次の桁で収まるとは限らず、さらに次の桁にまで持ち越す場合が出てきた。結果を以下に示す。桁位置を合わせるため空白を書き加え、さらに注釈も書き加えてある。

                                                 コンピュータでの数値表現の範囲
                                        1 =  1!
                                        2 =  2!
                                        6 =  3!
                                       24 =  4!
                                      120 =  5!   8ビット符号なし  
                                      720 =  6!
                                     5040 =  7!
                                    40320 =  8!  16ビット符号なし
                                   362880 =  9!   
                                  3628800 = 10!   
                                 39916800 = 11!   
                                479001600 = 12!  32ビット符号なし
                               6227020800 = 13!   
                              87178291200 = 14!   
                            1307674368000 = 15!   
                           20922789888000 = 16!   
                          355687428096000 = 17!   
                         6402373705728000 = 18!   
                       121645100408832000 = 19!   
                      2432902008176640000 = 20!  64ビット符号なし
                     51090942171709440000 = 21!   
                   1124000727777607680000 = 22!   
                  25852016738884976640000 = 23!   
                 620448401733239439360000 = 24!   
               15511210043330985984000000 = 25!   
              403291461126605635584000000 = 26!   
            10888869450418352160768000000 = 27!   
           304888344611713860501504000000 = 28!   
          8841761993739701954543616000000 = 29!   
        265252859812191058636308480000000 = 30!   
       8222838654177922817725562880000000 = 31!
     263130836933693530167218012160000000 = 32!
    8683317618811886495518194401280000000 = 33!
  295232799039604140847618609643520000000 = 34! 128ビット符号なし
10333147966386144929666651337523200000000 = 35!

もっと実用的なプログラムを書くなら、コンピュータの計算能力をより効率よく利用するものにすべきである。単純な改良としては基数を100とするか（対応して出力部での変換が必要になる）、または各種変数をより大きくして（例えば32ビット整数）さらに基数を大きく、例えば10,000にする。十進以外の基数から十進数出力への変換には多大な計算を要する。とはいうものの、コンピュータの本来扱える整数の限界に近い基数を使った方が有利である。通常の整数型であれば、その中身が6であろうと10000であろうと操作に要する時間は同じであり、大きい値を格納した方が配列全体としてより大きな数を表せる。コンピュータによっては、積をその桁の値と桁上り（キャリー）に自動的に分ける機能があり、例にある mod と div という操作を必要としない。例えば IBM 1130 では16ビット整数の乗算の積は32ビットで表され、これを2つの16ビット整数として扱うこともできる。その場合、大きな数の基数は65536となり、桁上りは上位16ビットで、桁は下位16ビットということになる。したがって、それらを分ける際に mod や div を必要としない。

プロセッサには存在するキャリーフラグなどを扱えない、などといった理由から、こういったたぐいのプログラミングには高水準言語は不適という面があり、機械語によるハックによって格段に高速な実装が実現できる場合もある。しかし、数値演算のコードはデバッグの難しさもあり、そのあたりのトレードオフもある。高水準言語のソースレベルでポータブルなライブラリが必要な場合、C言語では共用体、あるいは FORTRAN のEQUIVALENCE文や PL/1 の OVERLAY文などを使ったトリック的な書き方で、いずれも一般論としてはあまり薦められない手法とされているものではあるが、実装が可能な場合もある。

1つの桁の乗算において、一時変数は (radix - 1)² + carry を保持できなければならず、carry の最大値は (base - 1) である。IBM 1130 の場合、レジスタが32ビットなので、16ビットの演算の途中結果が16ビットで表せる範囲を超えても処理を続行できる。高級言語では、大きな数の配列の各要素（つまり桁）を16ビット符号無し整数とし、基数を65536としたとき、桁の乗算結果は 4,294,901,760 を超えないが、32ビット符号あり整数の上限は $2^{31}-1=2,147,483,647$ である。高級言語によっては、たとえハードウェアが32ビット符号無し整数の演算をサポートしていても、32ビット符号無し整数（上限は $2^{32}-1=4,294,967,295$ ）をデータ型として使えない場合がある。また、32ビット符号無し整数が使える場合でも、64ビット符号無し整数ではどうだろうか?

同様に、配列のインデックス用変数も16ビットや32ビットという制限がある。64ビット整数をインデックス変数に使えるとしても、今度はメモリ容量の限界という問題がある。さらに大きな数を表すには、適当な大きさのブロックに分け、インデックスとして（ブロック i、桁 j）のように2つの変数を使う方法や、インデックス自体も大きな数として桁ごとの配列で表す方法も考えられる。いずれにしても記憶装置容量の限界は避けられない。観測可能な宇宙に存在する全原子数は $10^{80}$ を超えないと予測されている。

歴史

1959年に発売された IBM 1620 は任意精度演算を実装した初期の例である。1620 は可変ワード長の十進コンピュータであり、メモリの許す限り任意の桁数の数値について計算が可能だった（浮動小数点の場合、仮数や指数の桁数には制限がある）。メモリを最大に搭載すると6万桁の数値を格納できるが、1620用FORTRANコンパイラは桁数を10桁に制限していた。1620はトランジスタを使っていたが、加算用の回路を持たず、メモリ上のテーブルを使って加算を実現していた。

IBM 初のビジネス用コンピュータ IBM 702 は511桁までの任意精度の数値の演算をハードウェアで実装していた。ソフトウェアでの初期の多倍長整数の実装としては、Maclisp が有名である。1980年代になると、VAX のオペレーティングシステムで数字の文字列を数値として計算する関数群が多倍長整数機能として用意されるようになった。また、 IBM では REXX が多倍長整数をサポートしていた。