数字列を多倍長整数のフォーマットに変換してくれる関数。

キー入力の際に必要となる。

<関数名>
  mpStr2Num ---- 数字列を多倍長整数に変換する

<形式>
  void mpStr2Num(int *num, char *str);

<引数>
  num (出力) 変換された多倍長整数
  str  (入力) 数字列

<関数値>
  なし

<注意事項>
  配列の各要素 ai(i は 1 以上)は1語を表し、1語で表し得る最大の整数は 9999 。語の長さは a0 の値で表す。すなわち、多倍長整数は

      anKn-1+an-1Kn-2+...+a2K+a1

   で表現する。ただし、K=10000、n=a0
  
用例

<関数本体>
  mpStr2Num.c

整数同士の比較。

<関数名>
  mpCmp ---- 多倍長整数同士の大小比較

<形式>
  int mpCmp(int *a, int *b);

<引数>
  a, b 多倍長整数

<関数値>
  aはbよりも大きい場合は正の値、小さい場合は負の値、等しい場合はゼロの値

<注意事項>
  配列の各要素 ai(i は 1 以上)は1語を表し、1語で表し得る最大の整数は 9999 。語の長さは a0 の値で表す。すなわち、多倍長整数は

      anKn-1+an-1Kn-2+...+a2K+a1

   で表現する。ただし、K=10000、n=a0
  
用例
  本関数以外に、キー入力に 関数 mpStr2Num()数字列を多倍長整数に)、画面表示に 関数 mpNum2Str()多倍長整数を数字列に)も使う。

<関数本体>
  mpCmp.c

<説明>
  語数の大きさをまず比べる。同じ語数ならば、上位語より順に比べていけばよい。

桁数の長い整数同士の割り算。

整数の世界から逃げ出す不思議な計算。実数の世界に通じるトンネルのようなものかな。計算結果としての有理数、無理数。言葉の遊びとも受け止められるかもしれないけど、大変深い洞察に基づく概念。

<関数名>
  mpDiv ---- 多倍長整数同士の除算

<形式>
  int mpDiv(int *q, int *r, int *a, int *b);

<引数>
  q   (出力) 多倍長整数同士の除算の商
  r   (出力) 多倍長整数同士の除算の余り
  a, b (入力) 多倍長整数 (a が被除数、b が除数)

<関数値>
  正常に計算できたときは 0。失敗したときは -1。除数が 0 のときは -2。

<注意事項>
  配列の各要素 ai(i は 1 以上)は1語を表し、1語で表し得る最大の整数は 9999 。語の長さは a0 の値で表す。すなわち、多倍長整数は

      anKn-1+an-1Kn-2+...+a2K+a1

   で表現する。ただし、K=10000、n=a0
  
用例
  本関数以外に、キー入力に 関数 mpStr2Num()数字列を多倍長整数に)、画面表示に 関数 mpNum2Str()多倍長整数を数字列に)も使う。

<関数本体>
  mpDiv.c

<説明>
  除算は四則演算の中に最もやっかいな計算である。ほかの演算と異なる点は、商の各桁の値の「見当をつける」という操作が含まれる点である。なるべく少なく手間で、適切な商の候補をみつけることが、プログラムの繁雑さを増してしまう。

多桁数の整数同士の掛け算についてだ。

<関数名>
  mpMul ---- 多倍長整数同士の乗算

<形式>
  void mpMul(int *ret, int *a, int *b);

<引数>
  ret (出力) 多倍長整数同士の乗算の結果(a * b)
  a, b (入力) 多倍長整数

<関数値>
  なし

<注意事項>
  配列の各要素 ai(i は 1 以上)は1語を表し、1語で表し得る最大の整数は 9999 。語の長さは a0 の値で表す。すなわち、多倍長整数は

      anKn-1+an-1Kn-2+...+a2K+a1

   で表現する。ただし、K=10000、n=a0
  
用例
  本関数以外に、キー入力に 関数 mpStr2Num()数字列を多倍長整数に)、画面表示に 関数 mpNum2Str()多倍長整数を数字列に)も使う。

<関数本体>
  mpMul.c

<説明>
  乗算は加減算ほど単純でなく、いろいろ工夫の余地がある。実用上または理論上の観点から興味深い方法がいくつも考えられている。

  乗算について考えよう。a、bをそれぞれm、n語とすると、積c=a*bはm+n 語以下となる。通常の筆算と同様の方法で、乗数bの1語ごとをaに掛け、それらの和は積cとなる。

  この方法では、bのn語をそれぞれaのm語に掛けるので、m*n回の繰り返しとなる。したがって、計算に必要な時間はO(mn)。

多桁数の整数同士の引き算。

<関数名>
  mpSub ---- 多倍長整数同士の減算

<形式>
  void mpSub(int *ret, int *a, int *b);

<引数>
  ret (出力) 多倍長整数同士の減算の結果 (a - b)
  a, b (入力) 多倍長整数 (a は b よりも大きいか等しい)

<関数値>
  なし

<注意事項>
   a は必ず b よりも大きいか等しいでなければならない。
  配列の各要素 ai(i は 1 以上)は1語を表し、1語で表し得る最大の整数は 9999 。語の長さは a0 の値で表す。すなわち、多倍長整数は

      anKn-1+an-1Kn-2+...+a2K+a1

   で表現する。ただし、K=10000、n=a0
  
用例
  本関数以外に、キー入力に 関数 mpStr2Num()数字列を多倍長整数に)、画面表示に 関数 mpNum2Str()多倍長整数を数字列に)も使う。

<関数本体>
  mpSub.c

<説明>
  10進数の計算と同じように、下語より引き算を実行する。上の語より「借り」が必要なときは、上の語に1を余分に引く。

どんなに桁数が長くても、整数同士の四則計算等を誤差なくやってくれる関数の紹介。

何百桁とか、長い整数を扱うのにいまのコンピュータは得意ではない。実数でも精度(いわゆる有効桁数)もそれほど高く保もたれていない。皆さんの普段使ってる電卓でもそう。液晶の表示部分は10桁とか、それぐらい。100桁も表示する電卓なんてないし、内部的にも有効桁数を沢山取って計算しているわけでもない。

しかし、われわれ普段習ってる数学では、数の有効桁数とか、使う範囲とかは考えない。理想的なものとしか扱わないから。無限大という数もふつうに使える。コンピュータでは絶対無理だけどね。

つまり、われわれ人間は数を概念的・記号的に思考するのは得意だが、ちょっとした計算、例えば、5,6桁の数同士の掛け算・割り算でもなかなか速くできないし、よく間違う。一方、コンピュータは計算が超得意だけど、記号とかになると苦手になることが多い。そう考えれば、人間とコンピュータとの合体、超高性能なコンピュータを腕時計や携帯のように、いつも持ち歩き、使うくことができれば、あるいは、人体に直に埋め込むことができれば、凄いことになりそう。

<関数名>
  mpAdd ---- 多倍長整数同士の加算

<形式>
  void mpAdd(int *ret, int *a, int *b);

<引数>
  ret (出力) 多倍長整数同士の加算の結果(a + b)
  a, b (入力) 多倍長整数

<関数値>
  なし

<注意事項>
  配列の各要素 ai(i は 1 以上)は1語を表し、1語で表し得る最大の整数は 9999。語の長さは a0の値で表す。すなわち、多倍長整数は

     anKn-1+an-1Kn-2+...+a2K+a1

  で表現する。ただし、K=10000、n=a0
  
用例
  本関数以外に、キー入力に 関数 mpStr2Num()数字列を多倍長整数に)、画面表示に 関数 mpNum2Str()多倍長整数を数字列に)も使う。

<関数本体>
  mpAdd.c

<説明>
  10進数の計算と同じように、下語より加えていく。和が10000 を超えたら桁上がりが発生し、上の語に1を余分に加える。

自然対数の底eの値を多桁数で算出する関数。

<関数名>
  e ---- 自然対数の底eの値を多桁数で算出する

<形式>
  void e(int *data, int keta);

<引数>
  data (出力) eの値が10進4桁ずつ入った整数の配列
  keta (入力) 小数点以下の桁数の指定

<関数値>
  なし

用例
  int data[1000];
  e(data, 1000);

  e の値 (小数点以下1万桁)

<関数本体>
  e.c

<説明>
  マクローリン展開

ex = 1 + x/1! + x2/2! + ... +xn/n! + ...

に x = 1 とおき、

     

e = 1 + 1/1! + 1/2! + 1/3! + ... + 1/n! + ...

の式を利用して計算。さらに、右辺の第2項目以降の計算では、計算回数の少ないつぎの式に変形して使う。

     (((...((1/n + 1) 1/(n-1) + 1) 1/(n-2) + 1)...) 1/3 + 1) 1/2 + 1

円周率 π を精度よく計算する関数。

<関数名>
  pi ---- 円周率 π の値を多桁数で算出する

<形式>
  int pi(int *data, int keta);

<引数>
  data (出力) 円周率の値が10進4桁ずつ入った整数の配列
  keta (入力) 小数点以下の桁数の指定

<関数値>
  正常に計算できたときは 0、エラーのときは -1。

用例
  int data[1000];
  pi(data, 1000);

  π の値 (小数点以下1万桁)

<関数本体>
  pi.c

<説明>
  πの値を求めるには、いろいろな式があるが、arctan(逆正接関数)を無限級数で表した式がよく使われる。式が単純な形で扱いやすく、ある程度の桁数の値を得るには適している。しかし欠点は、非常に長い桁数を得るには時間がかかりすぎることである。
  ここで使う式は、最もよく知られている、J. Machinの式
     π/4 = 4 arctan (1/5) - arctan (1/239)
  である。彼は1706年にこの公式でπを小数点以下100桁を求めた。同じ公式で、1949年にはENIACにより2037桁が約70時間で、1958年にはIBM704により1万桁が約100分で求められている。
  arctanをによる古典的な公式にほかに以下のものがある。
  ガウスの公式
     π/4 = 12 arctan (1/18) + 8 arctan (1/57) - 5 arctan (1/239)
  クリンジェンシェルナの公式
     π/4 = 8 arctan (1/10) - arctan (1/239) - 4 arctan (1/515)
  シュテルマーの公式
     π/4 = 6 arctan (1/8) + 2 arctan (1/57) + arctan (1/239)
  高野喜久雄の公式
     π/4 = 12 arctan (1/49) + 32 arctan (1/57) - 5 arctan (1/239) + 12 arctan (1/110443)
  柴田昭彦の公式
     π/4 = 17 arctan (1/22) + 3 arctan (1/172) - 2 arctan (1/682) - 7 arctan (1/5357)

  シュテルマーの公式に8の割り算を工夫すれば高速化期待できるし、高野の公式を利用して並列処理をすれば良い成績が期待できよう。

  arctanを使った公式でπを計算すると、桁数を2倍にするのに計算時間は4倍の増加になるが、以下の方法では2.1〜2.2倍で済む。

  ガウス・ルジャンドルの公式

     A = 1、B = (1/2)1/2、T = 1/4、X = 1

  として、AとBの差が必要とする精度より大きい間、次の計算を繰り返し実行する。

     Y = A、A = (A + B)/2、B = (B * Y)1/2
     T = T - X * (Y-A)2、X = 2 * X

  すると、πの値は(A + B)2/(4 * T) となる。ただし、A、B、T、Yの各値は、求めようとする精度以上で計算しておく必要がある。

1 から 99 の間の平方根の値を多桁数で算出する関数。多くの桁数の平方根が必要な時に利用するといいだろう。

平方根の算出は一般的に難しいけど、結果があっているかどうかの検証はものすごく簡単。2乗してみればすぐに判るからだ。円周率 π とは大違い。

では、関数の紹介。

<関数名>
  squareRoot ---- 1 から 99 の間の平方根の値を多桁数で算出する。

<形式>
  int squareRoot(int *data, int keta, int n);

<引数>
  data (出力) 平方根の値が10進4桁ずつ入った整数の配列
  keta (入力) 小数点以下の桁数の指定
  n   (入力) 計算しようとする平方根の2乗の値 (1 〜 99)

<関数値>
  正常に計算できたときは、計算結果となる data の長さ。エラーのときは 0。

用例
  2 の平方根を小数点以下1000桁で求める例
  int data[1000];
  squareRoot(data, 1000, 2);

  2 の平方根の値 (小数点以下1万桁)
  3 の平方根の値 (小数点以下1万桁)
  5 の平方根の値 (小数点以下1万桁)
  6 の平方根の値 (小数点以下1万桁)
  7 の平方根の値 (小数点以下1万桁)
  8 の平方根の値 (小数点以下1万桁)
  10 の平方根の値 (小数点以下1万桁)

<関数本体>
  squareRoot.c

<説明>
  必要とする桁数以上に、以下の漸化式

      pi+1 = pi2 - 2
      qi+1 = pi qi

  を計算し、最後に、pi+1 / qi+1 で平方根を求める。

  なお、初期値 p0、q0ペル方程式 p02 - M q02= 4 を満たす整数解。ただし、M は計算しようとする平方根の2乗の値。

  数学的帰納法で証明できるが、上の漸化式で得られた pi+1、qi+1 により、

        pi+1 / qi+1 = (M + 4 / (pi qi)2)1 / 2

  となる。したがって平方根の計算において、本アルゴリズムの方法ではニュートン法等ほかの方法よりも効率がよく、数回の計算で必要とする精度のものが得られる。

階乗 n!
     n! = 1 x 2 x 3...x n

をまともにと計算していくと、かなりのコンピュータでも、すぐにオーバフローしてしまう。ということで、本関数のようにそれなりの工夫が必要。

なお、n! の桁数は スターリングの近似式

n! = nne-n(2nπ)1/2

より、

   

n log10(n/e) + (log10(2nπ)/2+1

で求められる。

<関数名>
  factorial ---- 階乗 n! の値を多桁数で算出する

<形式>
  int factorial(int *data, int n);

<引数>
  data (出力) n!の値が10進4桁ずつ入った整数の配列
  n   (入力) 階乗 n の値

<関数値>
  配列dataに計算値の入った要素の数(配列の先頭から)

用例
  int data[1000];
  factorial(data, 1000);

  1000! の値

<関数本体>
  factorial.c