多倍長整数の四則演算 - ちりもつもればミルキーウェイ

はじめに

これは RSA完全理解 Advent Calendar の15日目の記事です。

RSAについてはひとしきり説明が終わったので、スクラッチで書くために巨大な数の四則演算について説明しようと思います。

四則演算さえできれば！いままでやったユークリッド互除法で逆元ポンするやつとかをシュッと実装できるぞ！盛り上がってきた

多倍長整数とは

ようはbigint的なものをイメージしてもらえばよいです。
以前のRSAの実装とかでもしれっと使ってましたね。

サイズの大きい数値は普通にint64とかで計算するとオーバーフローしちゃうので、なにか別の管理が必要です。

そこで、巨大な数値を一桁の数字のsliceとして管理します。これで四則演算を一桁の数値の演算で表現できれば良さそうです。

(なおこれは簡略化した実装であり、本来はシステムの単精度符号なし整数(64bit or 32bit)のuintで桁上りするような数列で管理したほうが無駄が少ないです)

で、巨大整数の演算っていうのは計算オーダーの影響をめちゃめちゃ受けやすいので、できる限り効率よく実装すべきです。(とくに乗算/除算については $\mathcal{O}(n \ log \ n)$ のアルゴリズムもあるらしいぞ。整数論的なフーリエ変換つかったりするらしくてめちゃくちゃむずそうだけど)

なんですが、いきなりちゃんと作っていろいろ説明とかすると大変だと思うので、まずはシンプルに筆算の考え方で愚直に実装したいと思います。筆算もあんまりナメられなくて、加算と減算は $\mathcal{O}(n)$ なので十分早い(なんなら定数倍とか除いたオーダーで比較するとこれ以上早いのないのでは)

まず符号なしの演算について実装し、つぎに符号付きの考慮をしたいと思います。

型定義とかコンストラクタとか

ひとまず型の準備をします。符号と一桁の数字のsliceをもつ型を定義します。
一桁の数字は一旦 uint8 とします。

type Int struct {
    neg bool
    abs digits
}

type digits []uint8

digits は最終的には一桁数字のscileになりますが、演算中に桁があふれる可能性があります。
今回の例では uint8 で足りるので一旦そのようにすすめます！

先程もちらっと触れましたが、本来は uint としつつ、数字1桁だけでなくその範囲の値すべて表現するようにしたほうが格段に効率が良いですが、簡単のため今回は uint8 の数字一桁のscileとして扱います

Int型のコンストラクタは以下のように書けます

func NewInt(x int64) *Int {
    neg := false
    if x < 0 {
        neg = true
        x = -x
    }
    abs := make(digits, 0)
    for i := x; i >= 1; {
        d := uint8(i % 10)
        if len(abs) == 0 {
            abs = append(abs, d)
        } else {
            abs = append([]uint8{d}, abs...)
        }
        i /= 10
    }
    return &Int{
        neg: neg,
        abs: abs,
    }
}

一桁ずつ10の剰余を求めてsliceに突っ込むだけですね。

注意としては、剰余とって見てる都合で下の桁から取ってるので、sliceにはどんどん先頭にpushするように書いてます。
これで実際に見える数字と同じ順番でsliceに入ります

あとは開発中にあると便利なので fmt.Stringer を実装してあげましょう。これはちゃんと書かずにサボります。標準パッケージばんざい

func (b *Int) String() string {
    var s string
    if b.neg {
        s += "-"
    }
    for _, v := range b.abs {
        s += fmt.Sprintf("%d", v)
    }
    return s
}

digitsについて、たとえば演算結果で繰り下がりが発生して最上位1桁が0となってしまう場合がありえるので、余計な0を消すための処理も追加します

func norm(abs digits) digits {
    i := 0
    for i < len(abs) && abs[i] == 0 {
        i++
    }
    return abs[i:]
}

あとdigits同士の比較もしたい場合があるので、そのような処理も追加しておきます

func cmp(x, y digits) int8 {
    // サイズが異なる場合サイズの大きいほうが絶対値が大きい
    m, n := len(x), len(y)
    if m != n {
        if m > n {
            return 1
        } else {
            return -1
        }
    }

    // サイズが同一の場合はより上位の桁の値が大きいほうが絶対値が大きい
    for i := 0; i < m; i++ {
        dx, dy := x[i], y[i]
        switch {
        case dx > dy:
            return 1
        case dx < dy:
            return -1
        }
    }
    return 0
}

よくある

左が大きければ1
左と右が等しければ0
右が大きければ-1

を返すような比較関数です。
digitsについてはこれくらいあれば大丈夫っすね！

符号なし加算

$789 + 789$ みたいな計算について、筆算で考えてみます

1の位は $9+9 = 18$ だから $10$ を次の桁に繰り上がりさせて $8$
10の位は $8+8+1=17$ だから $10$ を次の桁に繰り上がりさせて $7$
100の位は $7+7+1=15$ だから $10$ を次の桁に繰り上がりさせて $5$
1000の位は $0+0+1$ だから $1$
合計すると $1578$ となる！

これは各桁の位をそれぞれ下から足し合わせていけば簡単に実装できそうですね！
繰り上がりの管理だけ少し気をつければ特に難しいことはなさそうです。

繰り上がりを考慮して末尾の桁から計算するので、ループの向きもそれに揃えます。

func add(x, y digits) digits {
    m, n := len(x), len(y)
    var l int
    if m > n {
        l = m
    } else {
        l = n
    }
    abs := make(digits, l+1) // capは繰り上がり考慮で+1
    for i := l - 1; i >= 0; i-- {
        var dx, dy uint8
        if (m-1)-(l-1-i) >= 0 {
            dx = x[(m-1)-(l-1-i)]
        }
        if (n-1)-(l-1-i) >= 0 {
            dy = y[(n-1)-(l-1-i)]
        }
        // 繰り上がりを考慮した各桁の加算
        abs[i+1] = abs[i+1] + dx + dy
        if abs[i+1] >= 10 {
            // 結果が10以上だったら次の桁に繰り上がり
            abs[i+1] -= 10
            abs[i] = 1
        }
    }
    return norm(abs)
}

dx = x[(m-1)-(l-1-i)] となっているところが少しややこしいですが

m-1 はxの最後の要素のindex
l-1-i はループの最大値からループの現在値を引いた数 = 現在のステップ数

ということなので、slice xから現在のステップで参照すべき要素を末尾から抜いていってるってことです。(配列の要素を末尾から取り出すのに合わせて、ループの方向を大きい方から小さい方に回しているのですこし複雑になっている)

この処理の計算量は、足し合わせる数のうち桁数の大きい方の桁数に依存し、それぞれの桁について演算を行うので

$\mathcal{O}(n)$

といえそうです。まあかなり小さい計算量なのではないだろうか。
Goの mith/big の実装とかも読んでみましたが概ね同じことをしていたので、細かい書き方で定数倍は違うかもですがオーダーは似たようなもんだと思います。

符号なし減算

これも筆算で考えてみます。繰り下がりの $1234 - 987$ で考えてみます

1の位は $4-7$ ができないので、10の位から繰り下げて $14-7=7$ と整理 $7$
10の位は $(3-1) - 8$ ができないので、100の位から繰り下げて $12-8=4$ と整理して $4$
100の位は $(2-1)-9$ ができないので、1000の位から繰り下げて $11-9=2$ と整理して $2$
1000の位は $(1-1)-0$ で $0$

これも簡単に実装できそうですね。uint8だと負数が表現できなくて繰り下がりはちょっと頭を使いますが、まあそんなに複雑ではないです。

また、今回は符号なし減算について考えているので、たとえば $100-200$ のように引く数のほうが大きくなると取り扱いにこまります。なのでそのような入力が来たらpanicにしちゃいましょう。
いいのかそれはみたいな感じもありますが、あくまで符号を考えるときになんとかしてあげればよいです。

ちなみに math/big は非公開な符号なし減算処理について引く数のほうが多いときは同じくpanicにしてます。非公開処理なのでライブラリ内で使い方に気を配ればよいだろうという整理ですね。

func sub(x, y digits) digits {
    m, n := len(x), len(y)
    if m < n || cmp(x, y) == -1 {
        panic("underflow")
    }
    var l int
    if m > n {
        l = m
    } else {
        l = n
    }
    abs := make(digits, l)
    for i := l - 1; i >= 0; i-- {
        var dx, dy uint8
        if (m-1)-(l-1-i) >= 0 {
            dx = x[(m-1)-(l-1-i)]
        }
        if (n-1)-(l-1-i) >= 0 {
            dy = y[(n-1)-(l-1-i)]
        }

        // 繰り下がりを考慮した各桁の減算
        if dx < abs[i] || dx-abs[i] < dy {
            abs[i] = (dx - abs[i] + 10) - dy
            // 次の桁に繰り下がりのflagを立てておく
            // |x| >= |y| が保証されているので最上位桁の繰り下がりは考慮しなくてよい
            abs[i-1] = 1
        } else {
            abs[i] = (dx - abs[i]) - dy
        }
    }
    return norm(abs)
}

減算も加算と同じように桁数の大きい方の桁数に計算量が依存します。筆算のしくみを実装したときの減算の計算量は

$\mathcal{O}(n)$

です。加算とおなじく、書き方によって定数倍はちがうかもですがだいたいどの言語もこんなもんな印象があります

符号なし乗算

これも筆算のやりかたで考えます。

(はじめにいっておくと、乗算における筆算のアルゴリズムは比較的効率がよくなく実際のbigint的な実装ではまず使われません。ただ考え方がわかりやすいので、まずはここから実装します。ブラッシュアップは余裕があればやるぞ)

$789 + 789$ の例で考えます。筆算のやり方を整理すると

$9*9=81$ なので10の位に8, 1の位に1を置く。加算と違ってすべての桁について計算する必要があるのでこの時点では確定しない
$9*80=720$ なので100の位に7,10の位に2を置く
$9*700=6300$ なので1000の位に6, 100の位に3を置く
$80*9=720$ なので100の位に7, 10の位に2を置く
$80*80=6400$ なので1000の位に6, 100の位に4を置く
$80*700=56000$ なので10000の位に5, 1000の位に6を置く
$700*9=6300$ なので1000の位に6, 100の位に3を置く
$700*80=56000$ なので10000の位に5, 1000の位に6を置く
$700*700=490000$ なので100000の位に4, 10000の位に9を置く
すべて加算したものが答え

という段取りです。わかりやすいように表にすると

\	700	80	9
700	400000 + 90000	50000 + 6000	6000 + 300
80	50000 + 6000	6000 + 400	700 + 20
9	6000 + 300	700 + 20	80 + 1

となり622521が答えになります。

筆算の考え方を愚直に実装すると以下のようになります

func mul(x, y digits) digits {
    m, n := len(x), len(y)
    if m == 0 || n == 0 {
        return digits{}
    }
    abs := make(digits, m+n)
    for i := m - 1; i >= 0; i-- {
        for j := n - 1; j >= 0; j-- {
            dx := x[i]
            dy := y[j]
            // 繰り下がりを考慮した各桁の乗算
            abs[i+j+1] += dx * dy
            if abs[i+j+1] >= 10 {
                carry := abs[i+j+1] / 10
                abs[i+j+1] -= carry * 10
                abs[i+j] += carry
            }
        }
    }
    return norm(abs)
}

x, yの各桁について総当りで乗算して、その都度繰り上がり考慮する！といった感じです。

計算量について考えてみると、総当りで計算するので入力された数値両方の桁数に依存します。特に $n$ 桁 $\times$ $n$ 桁のとき

$\mathcal{O}(n^2)$

と言えそうです

符号なし除算

除算も乗算とおなじく筆算の手順でやるとすこし効率がわるいです。

なんですが、おなじく簡単のためにここから出発します。 $1000 \div 23$ の例で試します

$10000 \div 23$ について100の位から考えて $2300 \times 4 = 9200$ なので100の位に4, 計算途中のあまりに800を置く
$800 \div 23$ について10の位で考えて $230 \times 3 = 690$ なので10の位に3, 計算途中のあまりに110を置く
$110 \div 23$ について1の位で考えて $23 \times 4 = 92$ なので1の位に4, これ以上は桁を下げられないのであまりは8
$10000 \div 23$ の答えは $434 \dots 8$

この流れを実装してみましょう。除算は上の桁から計算していくので、ループの向きは今までと逆なので注意です。

除算は筆算界では一番複雑です。コメント丁寧に振りましたが少し難しいかも

func div(x, y digits) (quo digits, rem digits) {
    m, n := len(x), len(y)
    if n == 0 {
        panic("division by zero")
    }
    if cmp(x, y) < 0 {
        return digits{}, x
    }
    // 結果の桁数lを求める
    l := m - n // 少なくとも x > y*10^(m-n-1) で、そのとき商の桁数は m-n
    if cmp(x[:n-1], y) >= 0 {
        l += 1 // y*10^(m-n) したときにxのほうが大きければさらにその桁についても商が存在する
    }

    // 内部の乗算のコストをかるくするため、remは除算中の桁に関連するあまりだけ保持する
    // y*10^l のときの10の指数ぶんの桁数はそのステップでは値が変動しないので初期値から除いておく
    remPos := (m) - (l - 1)
    rem = x[:remPos]
    quo = make(digits, l)
    for i := 0; i < l; i++ {
        for j := 1; j < 10; j++ {
            // 現在のあまりを超えない範囲で最大の y * j (jは1-9) を求めて今の桁の商とする
            if cmp(mul(y, digits{uint8(j)}), rem) <= 0 {
                quo[i] = uint8(j)
            } else {
                break
            }
        }
        rem = sub(rem, mul(y, digits{quo[i]}))
        if i < l-1 {
            remPos += 1
            rem = append(rem, x[remPos-1])
        }
    }
    return quo, rem
}

ポイントは、商の桁数はすくなくとも 入力の桁数の差 以上であるということです。その上で上位の桁について割られる数のほうが大きければさらに1桁追加されます。

少しむずかしいのはこの部分くらいで、あとは愚直に行けるとこまで掛け算やってみて、いけたところがその桁の商だ！という脳筋戦法でやっています。

さて計算量についてですが、こやつは比較だったり乗算だったりを内部でやってるので定数倍がだいぶ大きいです。

入力された数の桁数に違いがあるほど、筆算で言う下の段におりる回数が増えてコストが増えます。また入力された桁数に違いがすくないとそれはそれで内部で乗算される桁数が増えてコストが大きいです。(定数倍とはいえ運悪いと10回乗算するからな...)

$n$ 桁と $n/2$ 桁の組み合わせあたりが最悪のケースで、めちゃめちゃ下の段におりなきゃいけなくてかつ内部でそこそこのサイズの乗算をこなす事になります。その時の計算量は

$\mathcal{O}(n^2)$

なようです。

符号あり演算！

符号なしの計算ができたら符号ありの対応はさほど難しくないです。単に加算や減算を絶対値の演算に置き換えてあげれば良いだけです。

加算は以下のように書けます

func Add(x, y *Int) *Int {
    neg := x.neg
    var abs digits
    if x.neg == y.neg {
        abs = add(x.abs, y.abs)
    } else {
        // xとyの正負が異なっていれば絶対値についての減算と言い換えることができる
        if cmp(x.abs, y.abs) >= 0 {
            abs = sub(x.abs, y.abs)
        } else {
            // (-x) + y == y - x となるので結果のsignを反転させれば成り立つ
            abs = sub(y.abs, x.abs)
            neg = !neg
        }
    }
    return &Int{
        neg: len(abs) > 0 && neg,
        abs: abs,
    }
}

ポイントは異なる符号が渡されたときに減算に言い換えているところです。 $x + (-y) = x - y$ だよなという話です。
あと0を符号なしの値として扱いたいのでneg判定のところに条件を入れてます。

こうすることで大した手間なく符号付きの演算に対応できました

つぎは減算についてです。

func Sub(x, y *Int) *Int {
    neg := x.neg
    var abs digits
    if x.neg != y.neg {
        // xとyの正負が異なれば絶対値についての加算と言い換えることができる
        abs = add(x.abs, y.abs)
    } else {
        if cmp(x.abs, y.abs) >= 0 {
            abs = sub(x.abs, y.abs)
        } else {
            // (-x) - (-y) == y - x となるので結果のsignを反転させれば成り立つ
            abs = sub(y.abs, x.abs)
            neg = !neg
        }
    }
    return &Int{
        neg: len(abs) > 0 && neg,
        abs: abs,
    }
}

こちらも考えは加算とほぼ同じで、減算は異なる符号が渡されたときに加算に言い換えることができます。 $x - (-y) = x + y$ だよなという話です。簡単だ！

次は乗算です。こっちのほうがもっと簡単です。

func Mul(x, y *Int) *Int {
    abs := mul(x.abs, y.abs)
    return &Int{
        neg: len(abs) > 0 && x.neg != y.neg,
        abs: abs,
    }
}

乗算は符号あってもなくても絶対値の演算としてはやることかわらないので、加算や減算とくらべて随分シンプルです。
また最終的な符号の判断も、乗算はそれぞれの符号がことなる場合のみ負となるので簡単です。

つぎは除算！これもほぼ乗算とおなじ！

func Div(x, y *Int) (quo *Int, rem *Int) {
    q, r := div(x.abs, y.abs)
    quo = &Int{
        neg: len(q) > 0 && x.neg != y.neg,
        abs: q,
    }
    rem = &Int{
        neg: len(r) > 0 && x.neg,
        abs: r,
    }
    return quo, rem
}

商とあまりの2つがあって、商は乗算と同じ条件で正負の判断ができます。

あまりは元の数と同じ符号だよねというのはまあそれはそうという感じなので、やっぱり簡単です。

これで巨大な整数の四則演算ができたぞ！！！

実際に使ってみた

というわけでオレオレbigint実装がこちらにあります！

github.com

これをつかってint64からあふれるような巨大な結果になる演算すっぞ！

package main

import (
    "fmt"
    "math"
    "math/big"

    mybig "github.com/convto/mycrypto/big"
)

func main() {
    a := big.NewInt(math.MaxInt64)
    b := big.NewInt(math.MaxInt64)
    stdMul := new(big.Int).Mul(a, b)
    fmt.Println(stdMul)
    m := mybig.NewInt(math.MaxInt64)
    n := mybig.NewInt(math.MaxInt64)
    myMul := mybig.Mul(m, n)
    fmt.Println(myMul)
}

$ go run main.go

85070591730234615847396907784232501249
85070591730234615847396907784232501249

わーいわーいうごいたぞー

playgroundも貼っときます！(module取得で時間がかかるので機嫌が悪いと時間切れします)

go.dev