任意精度演算

任意精度演算（にんいせいどえんざん）^[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 が多倍長整数をサポートしていた。

任意精度演算をサポートしているソフトウェア

アプリケーション

Mathematica などの数式処理システムは任意精度演算をサポートしている。Mathematica はGMPを利用。
PARI/GP — オープンソースの計算機代数ソフトウェアで、任意精度をサポート。
bc - Unix系システムに標準で搭載されている任意精度計算プログラム
Maxima — 数式処理システム。Common Lisp で実装されていたため、多倍長整数を受け継いでいる。
Ultra Fractal — フラクタル生成プログラム
Fractint — 各種フラクタルを生成・描画できる FLOSS ソフトウェア
Virtual Calc 2000 — 任意精度/基数の電卓（Windows用）
Calc — C言語スタイルの任意精度電卓（LGPL2）
SpeedCrunch — オープンソースでクロスプラットフォームの高精度電卓
TTCalc — オープンソースの任意精度電卓（TTMathを使用）
Big Online Calculator — 浮動小数点数の任意精度電卓。オンラインで利用。（TTMathを使用）
Full Precision Calculator — 有効桁数数万桁以上可能、計算可能桁数10^10^15程度までの高速な任意精度電卓。オンラインで利用。
高精度計算サイト『keisan』-フリー計算 — 有効桁数130桁まで、計算可能桁数10^10^7程度までの任意精度電卓。オンラインで利用。
多倍長電卓LM — 「数式入力型の多倍長電卓で初等関数も任意の桁数で計算可能分数計算も可能。絶対値が約10^100000000までの数値が扱えます。」とある。

ライブラリ

任意精度演算は多くの場合、専用ライブラリを呼び出すことで実装されている。そのライブラリがデータ型を定義し、数値を指定した精度で格納したり、計算を実行するサブルーチン群を提供している。

ライブラリによって数値の内部表現は異なる。整数のみを扱うライブラリもあるし、各種基数（十進、二進など）で浮動小数点数を格納するもある。分数（有理数）形式で数を格納するものもあれば、実数を全て表現できるとするものもある。

パッケージ/ライブラリ名	数値型	言語	ライセンス
apfloat	十進浮動小数点、整数、有理数、複素数	Java、C++	LGPLおよびフリーウェア
ARPREC and MPFUN	整数、二進浮動小数点数、複素二進浮動小数点数	C++、バインディングは C++ と FORTRAN	BSD
Base One Number Class	十進浮動小数点数	C++	有償
bbnum library	整数、浮動小数点数	アセンブリ言語、C++	改訂BSD
phpseclib	十進浮動小数点数	PHP	LGPL
BigDigits	自然数	C言語	フリーウェア[1]
Class Library for Numbers（CLN）	整数、有理数、複素数	C言語、C++	GPL
Computable Real Numbers	実数	Common Lisp
IMSL		C言語	商用
decNumber	十進数	C言語	ICU licence（MIT Licence）[2]
FMLIB	浮動小数点数	FORTRAN
GNU Multi-Precision Library（および MPFR（英語版））	整数、有理数、浮動小数点数	C言語、C++、他^[4]	LGPL
GNU Multi-Precision Library for .NET	整数	C#、.NET	LGPL
GNU Multi-Precision Library for Eiffel	整数、実数	Eiffel	LGPL
HugeCalc	整数	C++、アセンブリ言語	有償
IMath	整数、有理数	C言語	フリーウェア
IntX	整数	C#、.NET	改訂BSD
JScience LargeInteger	整数	Java
LibTomMath	整数	C言語、C++	パブリックドメイン
LiDIA	整数	C言語、C++	非商用利用については無料
MAPM	整数、十進浮動小数点数	C言語（バインディングは C++、Lua）	フリーウェア
MIRACL	整数、有理数	C言語、C++	非商用利用については無料
MPI	整数	C言語	LGPL
MPArith	整数、浮動小数点数、有理数	Borland Pascal、Delphi	zlib
mpmath	浮動小数点数、複素浮動小数点数	Python	改訂BSD
NTL	整数	C言語、C++	GPL
OpenSSL	整数	C言語	OpenSSL+SSLeay
TTMath library	整数、二進浮動小数点数	アセンブリ言語、C++	改訂BSD
vecLib.framework	整数	C言語	プロプライエタリ
W3b.Sine	十進浮動小数点数	C#、.NET	改訂BSD
Boost.Multiprecision	整数、浮動小数点	C++	Boost Software License

言語

組み込みまたは標準ライブラリの形式で任意精度演算をサポートするプログラミング言語を以下に挙げる。

Common Lisp — 多倍長整数、有理数、複素数をサポート
dc — POSIXの電卓ユーティリティ
Erlang — 組み込みの整数型が多倍長整数となっている。
Haskell — 組み込みの整数型が多倍長整数となっている。
Java — java.math.BigInteger クラス（整数）、java.math.BigDecimal クラス（十進整数）
.NET — System.Numerics.BigInteger 構造体（整数:4以降）
OCaml — Numライブラリで多倍長整数と有理数をサポート
Perl — bigintプラグマを指定することによりスコープ単位で多倍長整数指定が可能。bignumプラグマにより多倍長浮動小数点指定。それぞれMath::BitInt, Math::BitFloatモジュールで実装されている。
Python — 組み込みの int（第3.x版）、long（第2.x版）という整数型が多倍長整数。標準ライブラリには Decimal クラスもあり、桁数を指定可能。
REXX — Open Object Rexx や NetRexx も同様。
Ruby — 組み込みの Bignum という整数型が多倍長整数。標準ライブラリには BigDecimal クラスがあり、桁数を指定可能。
Scheme — 任意精度の整数と有理数をサポート（R⁵RS では推奨、R⁶RS では必須）
Scala — BigInt クラス.
Smalltalk — Squeak などの各種処理系でサポート。
JavaScript ー bigint 型（プリミティブな符号付き整数型）。通常の演算子を併用することで任意精度演算をサポートしている

脚注・出典

[脚注の使い方]

^ 英: arbitrary-precision arithmetic
^ 英: bignum arithmetic
^ R.K. Pathria 著、A Statistical Study of the Randomness amongst the first 10,000 Digits of Pi、1962年、Mathematics of Computation 誌、第16巻、N77-80、pp 188-97。「77」という数字の並びが1000桁のブロック内に28回出現することについて「Such an extreme pattern is dangerous even if diluted by one of its neighbouring blocks」と書いている。
^ GMPY

参考文献

Knuth, Donald, The Art of Computer Programming, ISBN 0-201-89684-2, Volume 2: Seminumerical Algorithms, Section 4.3.1: The Classical Algorithms
Implementing Multiple-Precision Arithmetic, Part 1
Nikolaevskaya et al. : "Programming with Multiple Precision", Springer-Verlag, ISBN 978-3-642-25672-1, e-ISBN 978-3-642-25673-8 (2012).
幸谷智紀 : 多倍長精度数値計算, 森北出版, ISBN 978-4-627-85491-8 (2019年11月)。
Fürer, M. : "Faster Integer Multiplication", SIAM J. Comput., Volume 39 Issue 3, 2009, pages 979–1005.
Harvey, D. and van der Hoeven, J. : "Integer Multiplication in Time O(nlogn)".
Karatsuba, A. : "The Complexity of Computations", Proceedings of the Steklov Institute of Mathematics, Volume 211, 1995, pages 169–183.
Schönhage, A. and Strassen, V. : "Schnelle Multiplikation grosser Zahlen", Computing, Volume 7, Issue 3–4, September 1971, pages 281–292.
Erica Klarreich : "Multiplication Hits the Speed Limit", Communications of the ACM, January 2020, Vol. 63 No. 1, Pages 11-13. doi=10.1145/3371387.

表話編歴データ型
ビット列	ビットトリットニブルオクテットバイトワードダブルワード（英）
数値	整数型符号付整数型十進型（英語版）有理数型（英語版）実数型複素数型固定小数点型浮動小数点型半精度単精度倍精度四倍精度八倍精度（英語版）拡張倍精度ミニフロート bfloat16 ブロック浮動小数点
ポインタ	物理アドレス型論理アドレス型（英語版）仮想アドレス型（英語版）参照型
テキスト	キャラクタ型ストリング型ヌル終端
複合	配列可変長配列連想配列構造体レコード共用体タグ共用体（英語版）タプルコンテナリストキュースタックセットツリー代数的データ型
その他	ブーリアン型 void型 null型列挙型再帰データ型トップ型（英語版）ボトム型関数の型（英語版）不透明型（英語版）シンボル型（英語版） Nullable型 Option型 Result型
関連項目	データ構造型システムプリミティブ型抽象型抽象データ型ボックス化動的束縛
カテゴリ

概要