ユークリッドの互除法の高速化

その数式、プログラムできますか？

ユークリッドの互除法の高速化の方法については様々存在するが、本書では以下のような方法で高速化していた。

github.com

１. 通常のやり方。aとbで大きい方から小さい方を減算していく。

int gcm0 (int a, int b)
{
    while (a != b) {
        if (b < a) { a = a - b; }
        else       { b = b - a; }
    }
    return a;
}

２. ちょっとif分の数を減らす。aからbを引けるだけ引いていって、引けなくなったら変数の位置をスワップすることでコードの再利用。

#define SWAP(a, b) { int tmp = a; a = b; b = tmp; }

int gcm1 (int a, int b)
{
    while (a != b) {
        while (b < a) { a = a - b; }
        SWAP(a, b);
    }
    return a;
}

３. 引き算のアルゴリズムを改良する。単純に考えれば、

        while (b < a) { a = a - b; }

の部分を表現すると、

int segment_remainder (int a, int b)
{
    while (b < a) { a = a - b; }
    return a;
}

なのだが、ここで効率化のために以下の定理を利用している。

r=segment_remainder(a,2b)である場合、以下の式が成り立つ。
* r<=b の場合 segment_remainder(a,b) = r
* r> b の場合 segment_remainder(a,b) = r-b

これをどのようにして活用しているかと言うと、

int fast_segment_remainder (int a, int b)
{
    if (a <= b) { return a; }                 // aとbの剰余だから、a<=bだったら、そりゃ剰余はaである
    if (a - b <= b) { return a - b; }         // aからbを一回引くだけでbよりも小さくなるのだったら、a-bが剰余である
    a = fast_segment_remainder (a, b + b);    // ↑よりもさらに引き算をする必要があれば、2bだけ引き算してみる
    if (a <= b) return a;                     // 上の公式の「r<=bの場合」。ただしrではなくaになっていることに注意！
    return a - b;                             // 上の公式の「r> bの場合」。ただしrではなくaになっていることに注意！
}

int fast_segment_gcm (int a, int b)
{
    while (a != b) {
        a = fast_segment_remainder (a, b);
        SWAP(a, b);
    }
    return a;
}

という訳で、それぞれa=1071,a=1029で実行してみたときの命令数は以下。

• gcm0 : 189命令
• gcm1 : 109命令
• fast_segment_gcm : 72命令

あれ、gcm0とgcm1でそんなに命令数変わる？と思ったらgcm1の方がコアになるループが小さいので速いのか。納得。

800001bc <gcm0>:
800001bc:       1085000a        beq     a0,a1,800001e8 <gcm0+0x2c>
800001c0:       00000000        nop
800001c4:       00851023        subu    v0,a0,a1
800001c8:       00a4182a        slt     v1,a1,a0
800001cc:       00a43023        subu    a2,a1,a0
800001d0:       0083100a        movz    v0,a0,v1
800001d4:       00c3280a        movz    a1,a2,v1
800001d8:       14a2fffa        bne     a1,v0,800001c4 <gcm0+0x8>
800001dc:       00402021        move    a0,v0
800001e0:       03e00008        jr      ra
800001e4:       00000000        nop
800001e8:       03e00008        jr      ra
800001ec:       00801021        move    v0,a0

800001f0 <gcm1>:
800001f0:       1085000f        beq     a0,a1,80000230 <gcm1+0x40>
800001f4:       00000000        nop
800001f8:       00a4182a        slt     v1,a1,a0
800001fc:       10600006        beqz    v1,80000218 <gcm1+0x28>
80000200:       00801021        move    v0,a0
80000204:       00451023        subu    v0,v0,a1
80000208:       00a2182a        slt     v1,a1,v0
8000020c:       1460fffe        bnez    v1,80000208 <gcm1+0x18>
80000210:       00451023        subu    v0,v0,a1
80000214:       00451021        addu    v0,v0,a1
80000218:       10a20003        beq     a1,v0,80000228 <gcm1+0x38>
8000021c:       00a02021        move    a0,a1
80000220:       1000fff5        b       800001f8 <gcm1+0x8>
80000224:       00402821        move    a1,v0
80000228:       03e00008        jr      ra
8000022c:       00000000        nop
80000230:       03e00008        jr      ra
80000234:       00801021        move    v0,a0