見出し画像

Excelと多倍長整数の四則演算の2

Excelで、とても大きな整数の計算をしてみよう!の2回目です。

16桁以上の数値だと表記がバグるExcelで、どのようにして大きい数字を扱うのかは、前回説明しています。ので、そちらをご覧ください。

今回は「掛け算」です。
それではいってみましょう!


第3章 乗算

掛け算の方法

前回も参照しましたQiitaの記事では、掛け算の方法として4つ説明しています。

  1. 筆算

  2. カラツバ法

  3. Toom-Cook法

  4. FFT(高速フーリエ変換)

「筆算」は小学校でならう大きい数の掛け算の方法です。
もう忘れたという方もいるかもしれませんので、復習しておきましょう。

なつかしの筆算

下の数を一桁ずつばらして、上の数と掛けながら、同時に繰り上がり処理も行って、最後は足し算をする。
これをどうやって、プログラムに落とし込みましょうか。
繰り上がりの処理は前回作りましたので、繰り上がりしないで掛け算してみましょう。

繰り上がり処理を省いた筆算

これなら配列の掛け算をした後に、斜め足し算をすればできそうです。
完成したコードがこちら

FuncName : _AnMulWc

Discription : 配列数の乗算。筆算。

Argments : ax, ay

Function definition :
=LET(
    _x, ax,
    _y, ay,
    _lx, COLUMNS(_x),
    _ly, COLUMNS(_y),
    _lz, _lx + _ly - 1,
    _mx, _x * TRANSPOSE(_y),
    _z, IF(
        _ly <= 1,
        _mx,
        REDUCE(
            EXPAND(CHOOSEROWS(_mx, 1), , _lz, 0),
            SEQUENCE(_ly - 1, , 2),
            LAMBDA(_a, _s,
                LET(
                    _sr, EXPAND(HSTACK(SEQUENCE(, _s - 1, 0, 0), CHOOSEROWS(_mx, _s)), , _lz, 0),
                    _a + _sr
                )
            )
        )
    ),
    _AnCarryAndFix(_z)
)

簡単に説明します。

  • 1~4行目は下準備です。

  • 5行目で掛け算をしています。「 _y 」をTRANSPOSE関数で縦に並べ替えて、ただ掛けているだけです。これで各桁を掛けた数字がマトリックス状にならびます。

配列の掛け算
  • 6行目のIF関数の中で、斜め足し算をしています。

  • 最後に_AnCarryAndFix関数で繰り上がり処理を行っておしまいです。

斜め足し算の良い方法を思いつかず、結局はREDUME関数でぐるぐる回す形になりました。
実行してみた結果たこちら

10桁×10桁の掛け算(筆算)

検算もあっているようです。
無事完成しました。
速度はどうでしょうか?

CnMulWc関数の速度チェック

500桁×500桁で1.76秒。平均2秒ぐらいかかっているようです。
やはりREDUCE関数でぐるぐるしているところが重そうですね。


追記 2025/2/1

斜め足し算の良い方法を思いついたので、追記します。
まずはコードです。

Function Name : _AnMulWc

Discription : 配列数の乗算。筆算方式

Argments : ax, ay

Function definition :
=LET(
    _x, ax,
    _y, ay,
    _lx, COLUMNS(_x),
    _ly, COLUMNS(_y),
    _lz, _lx + _ly,
    _m1, _x * TRANSPOSE(_y),
    _m2, EXPAND(_m1, , _lz, 0),
    _m3, TOROW(_m2),
    _m4, WRAPROWS(_m3, _lz - 1, 0),
    _z, BYCOL(_m4, LAMBDA(_m, SUM(_m))),
    _AnCarryAndFix(_z)
)

7~9行目が斜め足し算です。

  • 6行目(_m1)で掛け算をして、マトリックスを作ります。

  • 7行目(_m2)は、列を_lx + _ly に増やし、追加した分は0でパディングしています。

  • 8行目(_m3)は、マトリックスを1行配列に変換しています。

  • 9行目(_m4)は、_m3を_lx + _ly - 1列で折返し、再びマトリックスにしています。(最終行の足りない分は0パディング)

  • 最後にBYCOL関数を使って、各列ごとにSUMしています。

こちらの新関数を使った速度チェックがこちら。

CnMulWc関数の速度チェック(リニューアル)

2秒ぐらいかかっていた500桁×500桁が、0.13秒で終わりました。早い!


カラツバ法

では、カラツバ法を試してみましょう。
詳細は参照の記事を見ていただくとして、仕組みを簡単に説明すると、

  1. 例えば、$${A = 12345678, B = 98765432}$$とする。

  2. $${A = a_1×10000 + a_2, B = b_1×10000 + b_2 }$$とする

    1. つまり、$${a_1 = 1234, a_2=5678, b_1 = 9876, b_2 = 5432}$$

    2. また $${D = 10000}$$ とすれば、$${A = a_1×D+a_2, B=b_1×D+b_2}$$

  3. $${A×B = (a_1×D+a_2) × (b_1×D+b_2) = a_1×b_1×D^2 + (a_1×b_2 + a_2×b_1)×D + a_2×b_2}$$

  4. ここで、$${c_1 = a_1×b_1, c_2 = a_1×b_2 + a_2×b_1, c_3 = a_2×b_2}$$と置くと、

  5. $${c_2 = (a_1 + a_2) × (b_1 + b_2) - a_1×b_1 - a_2×b_2 =  (a_1 + a_2) × (b_1 + b_2) - c_1 - c_3}$$

  6. $${c_4 =  (a_1 + a_2) × (b_1 + b_2)}$$ とおけば、

  7. $${A×B = c_1×D^2 + (c_4 - c_1 - c_3)×D + c_3}$$

8桁×8桁の掛け算は、4桁の数字にばらしたところで、4桁×4桁の掛け算を4回行う必要がある(3番目の式)けれど、ちょっとした工夫で、3回の掛け算に落とし込んで計算量を減らそうという作戦ですね。すばらしい。
これをコードにしたのがこちら

FuncName : _AnMulKm

Discription : 配列数の乗算。カラツバ法を使用。

Argments : ax, ay

Function definition :
=LET(
    _x, ax,
    _y, ay,
    _lx, COLUMNS(_x),
    _ly, COLUMNS(_y),
    _z, IF(
        _lx + _ly < 16,
        _CnToAn(VALUE(_AnToCn(_x)) * VALUE(_AnToCn(_y))),
        LET(
            _sp, ROUNDUP(MAX(_lx, _ly) / 2, 0),
            _x1, _AnCarryAndFix(IF(_lx > _sp, TAKE(_x, , _sp), _x)),
            _x2, _AnCarryAndFix(IF(_lx > _sp, DROP(_x, , _sp), 0)),
            _x3, _AnAdd(_x1, _x2),
            _y1, _AnCarryAndFix(IF(_ly > _sp, TAKE(_y, , _sp), _y)),
            _y2, _AnCarryAndFix(IF(_ly > _sp, DROP(_y, , _sp), 0)),
            _y3, _AnAdd(_y1, _y2),
            _z1, _AnMulKm(_x1, _y1),
            _z2, _AnMulKm(_x2, _y2),
            _z3, _AnMulKm(_x3, _y3),
            _z4, _AnAdd(_z1, _z2),
            _z5, _AnAdd(_z3, -1 * _z4),
            _z6, _AnCarryAndFix(HSTACK(SEQUENCE(, _sp, 0, 0), _z5)),
            _z7, _AnCarryAndFix(HSTACK(SEQUENCE(, _sp * 2, 0, 0), _z2)),
            _z8, _AnAdd(_z1, _z6),
            _AnAdd(_z7, _z8)
        )
    ),
    _AnCarryAndFix(_z)
)

簡単に説明します。
まず、再帰を使っています。
どんどん分割していけば計算が楽になるからですね。
どこまで分割するのかを、LET関数内の6行目で判別しています。

  • 1~4行目は下準備

  • 5行目のIF関数で、再帰のおしりを決めています。

  • 6行目の「 _lx + _ly < 16 」は、桁数が足して16以下なら再帰終了という意味です。なぜ16かというと15桁までならExcelくんは普通に計算してくれるからです。

  • 7行目で配列数⇒文字列数⇒数値に戻して、普通に掛け算し、配列数に戻しています。

  • 9行目は、「_x」と「_y」の内、大きい方の半分の値をとっています。割り切れない場合は繰り上げます。

  • 10~12行目で「_x」を分割及び足し算して3つの子データを作成

  • 13~15行目は同じことを「_y」で行います。

  • 16~18行目でそれぞれを掛け算した値を作成し、

  • 19~24行目で、それらを足したり引いたりシフトしたりして答えを出しています。

  • 27行目で繰り上がり処理をして終了です。

では、実行結果を見てみましょう。

10桁×10桁の掛け算(カラツバ法)

問題ないようですね。
では速度をみてみましょう。

CnMulKm速度チェック

ふむ、500桁×500桁で0.5秒ぐらい、筆算の3~4倍は早いようです。
しかし、これ、落とし穴がありました。

150桁×150桁での速度比較

なんと、150桁ぐらいだったら、筆算の方が早いようです。
どうやら、足し算も150回を超えたらバカならなくなるということでしょうか。
斜め算の高速化をなんとかしたいですね。だれか教えてください。

Toom-Cook法

Toom-Cook法は、カラツバ法の発展形です。
が、ここではやりません。
なぜかというと、「割り算」を使うからです。
割り算が完成したあとに試してみたいと思います。

FFT(高速フーリエ変換)

さて、もっとも早いのがFFTを使った手法だそうです。
しかし、こちらは難しい!
試しにvunuさんのこちらの記事

の、基数4FFTを使わせてもらって、試してみたのですが、誤差が出るんですよね・・・。
最新の記事、

で、誤差なくできているようですので、こちらを勉強しなおして、再度チャレンジしたいと思います。

今回の勉強結果はこちら

デバッグよろしくおねがいします。

(追記:2025/1/26)
vunuさんの助言もあって、なんとか成功したので、続きをかきます。

まずは、FFT(高速フーリエ変換)を使った掛け算の仕組みについて。
いわゆる「畳み込み」というものなんですね。
Youtubeの「3 Blue 1 Brown Japan」の動画がためになったので知らない方はこちらをどうぞ。

上の動画内のキャプチャー画像を下に貼ります。

3Blue1BlownJapan. “畳み込みの仕組み | Convolution”. YouTube, 26 Jan. 2024

簡単に解説しますと、

  1. 配列「a」「b」があり、それぞれ長さが「n」「m」となっている

  2. それぞれをFFT(高速フーリエ変換)する

    1. この際、2つとも長さを「n + m」にそろえる

    2. 足りない分は0で埋めてからFFTすれば出力値の長さをかえられる

  3. 出力値「$${\hat{a}}$$」「$${\hat{b}}$$」を掛け算する

  4. 「$${\hat{a}・\hat{b}}$$」をIFFT(逆高速フーリエ変換)する

となります。
さて、理屈はさっぱりだけれども、やり方はわかった!
でも高速フーリエ変換のコードなんてどう書けばいいの?と思っていたらタイミングよくvunuさんが作っていました。感謝!!
というわけで、vunuさんのコードをコピーして作ったのがこちら

Function Name : _AnMulFft

Discription : 配列数の乗算。FFTArgments : ax, ay

Function definition :
=LET(
    _lx, COLUMNS(ax),
    _ly, COLUMNS(ay),
    _lz, 2 ^ ROUNDUP(LOG(_lx + _ly - 1, 2), 0),
    _x, TRANSPOSE(EXPAND(ax, , _lz, 0)),
    _y, TRANSPOSE(EXPAND(ay, , _lz, 0)),
    _fx, _FftFR(_x),
    _fy, _FftFR(_y),
    _fz, REDUCE(
        IMPRODUCT(INDEX(_fx, 1), INDEX(_fy, 1)),
        SEQUENCE(_lz - 1, , 2),
        LAMBDA(_a, _s, VSTACK(_a, IMPRODUCT(INDEX(_fx, _s), INDEX(_fy, _s))))
    ),
    _ifz, _IFftFR(_fz),
    _z, TRANSPOSE(ROUND(_ifz, 0)),
    _AnCarryAndFix(_z)
)

簡単に説明します。

  • 1~5行目は下準備

    • 1~3行目で長さを決めます。

    • 4,5行目は足りない分をEXPAND関数にて0で埋めて、TRANSPOSE関数で縦横を逆にしています。(vunuさんのコードが縦配列仕様なので

  • 6,7行目でFFTしています。

  • 8~12行は掛け算です。

    • IMPRODUCT関数(複素数の掛け算)にはスピルが使えなかったので、REDUCE関数でループさせてますが、もっと楽な方法ある?

  • 13行目で逆高速フーリエ変換

  • 14行目は誤差補正しながら、縦横を逆にしています。

  • 15行目で繰上げ処理をしておしまいです。

_FftFR関数と、_IFftFR関数の中身の話はここではしません。
最下部にファイルを置いておきますので、気になる方はダウンロードしてみてください。
では、出力結果です。

20桁×20桁(FFT)

無事、誤差なく表示できました。めでたい!
では、最速とうたわれるFFTの速さを見てみましょう。

1000桁×1000桁で速度比較

1000桁程度でもカラツバ法と比べて倍ぐらい早いようです。
さすがですね。

さて、では掛け算はいったん終了です。
次は割り算!
今回の勉強結果を置いておきます。

デバッグ宜しくお願いします。


追記(2025/2/1)

筆算と同じで、BYROW関数を使えば、同じ行同士の複素数掛け算ができることが分かったので修正します。

Function Name : _AnMulFft

Discription : 配列数の乗算。FFTArgments : ax, ay

Function definition :
=LET(
    _lx, COLUMNS(ax),
    _ly, COLUMNS(ay),
    _lz, 2 ^ MAX(ROUNDUP(LOG(_lx + _ly - 1, 2), 0), 3),
    _x, TRANSPOSE(EXPAND(ax, , _lz, 0)),
    _y, TRANSPOSE(EXPAND(ay, , _lz, 0)),
    _fx, _FftFR(_x),
    _fy, _FftFR(_y),
    _fz, BYROW(HSTACK(_fx, _fy), LAMBDA(_m, IMPRODUCT(_m))),
    _ifz, _IFftFR(_fz),
    _z, TRANSPOSE(ROUND(_ifz, 0)),
    _AnCarryAndFix(_z)
)

_fz の計算が1行で済みました。
では前回のコードと速度を比べてみましょう。

改良したFft版の速度!

やはり早いですね。
新しい筆算版よりも早い。
ではカラツバ法も含めて1000桁計算してみます。

筆算 vs カラツバ法 vs FFT

カラツバ法よりも圧倒的にFFTが早くなりました。
ところが、こんどは筆算がエラーになりました。
まぁ、これは配列オーバーフローですね。
新方式だと(1000 + 1000 ) × 1000 = 2,000,000 配列が必要なわけですが、Excelの配列上限が 2^20 = 1,048,576 なのであふれちゃったと。
でも筆算はそこそこ早いので、少ない数字でFFtと比べてみました。

筆算とFFTの速度比較

100桁ぐらいまでは筆算が勝ったりしますが、そこを超えるとFFTの圧勝になるようです。
100桁以下でもこのぐらいの速度差ならFFT版をメインにしてもよさそうですね。
では、修正版のファイルを置いておきます。

デバッグよろしくお願いします。

いいなと思ったら応援しよう!