任意精度演算

任意精度演算とは、数値の精度を必要ならいくらでも伸ばしたりできるような演算システムによる演算である。

1. 概要
-----

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

多倍長演算とも言うが、プログラミング言語によっては、多倍長整数  の名前が bignum であることもある。

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

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

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

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

3. 用途
-----

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

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

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

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

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

4. 実装上の問題
---------

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

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

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

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

加算や減算の最も単純なアルゴリズムは、単純に桁ごとに順次加算・減算していき、繰り上げや繰り下げを必要に応じて行うものである。この計算量は である。

乗算の場合、最も単純なアルゴリズムは筆算の計算方法をそのまま採用する手法で の計算量である。計算量が の乗算アルゴリズムとして、高速フーリエ変換に基づくショーンハーゲ・ストラッセン法がある。また、カラツバ法では の計算量である。 があまり大きくない場合、カラツバ法の方が性能がよい。

5. 例
----

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

大きな階乗の正確な値を求めたい場合、特別なソフトウェアが必要となる。以下の擬似コードは、1、1*2、*3、*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  ; begin columns.length := 1 ; columns.values [columns.length] := 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 ;   end end ;

procedure writeColumns  ; begin if 36 < columns.radix then raise OutOfRangeError ; const DIGITS : array of character = ['0'..'9','A'..'Z'] ;  for var ci := columns.length downto 1 do write  ;  end ;

begin for var n := 1 to 35 do begin var columns : COLUMNS_TYPE ; columns.radix := 10 ;  factorial  ; writeColumns  ; writeln  end end.<WtXmlEndTag />

この例で、詳細を見てみよう。最も重要なのは、大きな数の表現方法の選択である。この場合、階乗の計算なので整数だけでよく、固定小数点方式で事足りる。基数の冪は0から始まって大きくなっていくので、配列の先頭が0で、後ろの方ほど基数の冪が大きくなるという表現が便利である。配列のインデックスの始点や大きさの指定方法は言語によって異なるが、一般に計算の必要条件には関係しないので、ここでは単純化するために配列の始点を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とするか、または各種変数をより大きくしてさらに基数を大きく、例えば10,000にする。十進以外の基数から十進数出力への変換には多大な計算を要する。とはいうものの、コンピュータの本来扱える整数の限界に近い基数を使った方が有利である。通常の整数型であれば、その中身が6であろうと10000であろうと操作に要する時間は同じであり、大きい値を格納した方が配列全体としてより大きな数を表せる。コンピュータによっては、積をその桁の値と桁上りに自動的に分ける機能があり、例にある mod と div という操作を必要としない。例えば IBM 1130 では16ビット整数の乗算の積は32ビットで表され、これを2つの16ビット整数として扱うこともできる。その場合、大きな数の基数は65536となり、桁上りは上位16ビットで、桁は下位16ビットということになる。したがって、それらを分ける際に mod や div を必要としない。

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

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

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

6. 歴史
-----

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

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

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

7.1. アプリケーション
-------------

en などの数式処理システムは任意精度演算をサポートしている。en はGMPを利用。
en — オープンソースの計算機代数ソフトウェアで、任意精度をサポート。
bc - Unix系システムに標準で搭載されている任意精度計算プログラム
en — 数式処理システム。en で実装されていたため、多倍長整数を受け継いでいる。
en — フラクタル生成プログラム
en — 各種フラクタルを生成・描画できる en ソフトウェア
[1] — 任意精度/基数の電卓
[2] — C言語スタイルの任意精度電卓
[3] — オープンソースでクロスプラットフォームの高精度電卓
[4] — オープンソースの任意精度電卓
[5] — 浮動小数点数の任意精度電卓。オンラインで利用。
多倍長電卓LM — 「数式入力型の多倍長電卓で初等関数も任意の桁数で計算可能 分数計算も可能。絶対値が約10^100000000までの数値が扱えます。」とある。

7.2. ライブラリ
----------

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

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



7.3. 言語
-------

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

en — 多倍長整数、有理数、複素数をサポート
dc — POSIXの電卓ユーティリティ
da — 組み込みの整数型が多倍長整数となっている。
en — 組み込みの整数型が多倍長整数となっている。
en — java.math.BigInteger クラス、java.math.BigDecimal クラス
en — System.Numerics.BigInteger 構造体
fr — [6]ライブラリで多倍長整数と有理数をサポート
en — bigintプラグマを指定することによりスコープ単位で多倍長整数指定が可能。bignumプラグマにより多倍長浮動小数点指定。それぞれMath::BitInt, Math::BitFloatモジュールで実装されている。
en — 組み込みの int、longという整数型が多倍長整数。標準ライブラリには Decimal クラスもあり、桁数を指定可能。
en — en や en も同様。
en — 組み込みの Bignum という整数型が多倍長整数。標準ライブラリには BigDecimal クラスがあり、桁数を指定可能。
en — 任意精度の整数と有理数をサポート
en — BigInt クラス.
en — en などの各種処理系でサポート。
JavaScript ー bigint 型。通常の演算子を併用することで任意精度演算をサポートしている

8. 脚注・出典
--------

9. 参考文献
-------

Knuth, Donald, //The Art of Computer Programming,// ISBN 0-201-89684-2, Volume 2: Seminumerical Algorithms, Section 4.3.1: The Classical Algorithms
[7]
Nikolaevskaya et al. : "Programming with Multiple Precision", Springer-Verlag, ISBN 978-3-642-25672-1, e-ISBN 978-3-642-25673-8 .
幸谷智紀 : 多倍長精度数値計算, 森北出版, ISBN 978-4-627-85491-8 。
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".
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.

10. 関連項目
--------

除算