大数除去アルゴリズム

8198 ワード
C言語テクニック
概要
実際のプロジェクトでは,同僚がアルゴリズムを移植する際に64ビット整数の除算演算を行うことに遭遇した.探してみると、Linuxカーネルにはこの演算をサポートする関数do_div()があり、この関数はLinux/arch/arm/include/asm/div 64にある.hファイルで実装する.具体的な実現方法がよく分からないので、自分で大数を除いたアルゴリズムを書くことができるかどうか考えました.以下はアルゴリズムの内容ですが、不足点があれば、ご指摘ください.
注意:次の式およびコードでは、m:被除数n:除数q(quotient):商r(remainder):残高
以下、いずれも100/7を例として説明する.最終的なコード実装は文章の最後に置かれる.
継続的減算
m/nの商がq,残数がrであると仮定すると,この式:m=n*q+rを得ることができ,変換してm−n*q=r(r例えば100/7です.100 - 7*1 = 93 > 7 100 - 7*2 = 93 - 7 = 86 > 7 100 - 7*3 = 86 - 7 = 79 > 7 100 - 7*4 = 79 - 7 = 72 > 7 …… …… 100 - 7*14 = 9 - 7 = 2 < 7
これはrdiv1()関数である.
このアルゴリズムは2つの大数を除いた結果を算出できるに違いないが、mが非常に大きく、nが非常に小さい場合、例えばm=0 xffffffffffffffff;n=2で、この時悲劇になって、全部で0 x 7 fffffffffffffの減算をします.だからこのアルゴリズムは必ずしも効率的ではない.
バイナリ特性による継続減算
前は7*1、7*2、7*3、...、7*14を試し続けましたが、このような効率は比較的低いので、乗算を変えてみることはできません.すなわち、7*1、7*2、7*4、7*8、さらに乗せれば7*16=112>100なので、このときは100-7*8=44で次のループに続きます.100=7*8+7*6+2です.このアルゴリズムの実装はdiv2()関数である.
にぶんほう
除算操作を乗算操作に変換します.
7*100=700>100 7*50=350>100 7*25=175>100 7*12=84<100、このときq=12 7*(12+6)=126>100 7*(12+3)=105>100 7*(12+1)=91<100、このときq+=1、すなわちq=13 7*(13+1)=98<100、このときq+=1、すなわちq=14 7*(14+1)=105>100
2分が1のとき、最後の結果がmより大きく、最後から2回目の結果がmより小さいことを満たす限り、ループは終了する.
このアルゴリズムは実際には大数の除算とは言えない.乗算操作が使われているからだ.mとnが同時に大きな数であれば、例えばmとnが64ビット整数であれば、両者を乗算した結果は64ビットを超えるに違いない.したがって,このアルゴリズムは小さい数で除去するのに適している.ここで挙げたのは,演算の効率が高いことを重視しているからである.このアルゴリズムの実装はdiv3()関数である.
インプリメンテーションコード
/*
     div1()/div2()/div3()         ：
m       :    
n       :   
*remain :   

        。

 ：            n   0       。
*/

#include 

typedef unsigned long uint64_t;
typedef long int64_t;

uint64_t div1(uint64_t m, uint64_t n, uint64_t *remain)
{
    uint64_t quot = 0;
    uint64_t value = 0;

    if(m < n)
    {
        quot = 0;
        *remain = m;
    }
    else if(m == n)
    {
        quot = 1;
        *remain = 0;
    }
    else
    {
        value = m-n;
        quot = 1;
        while(value >= n)
        {
            quot++;
            m = value;
            value = m-n;
        }
    }

    *remain = value;
    return quot;
}


uint64_t div2(uint64_t m, uint64_t n, uint64_t *remain)
{
    uint64_t quot = 0;
    uint64_t value = m;
    uint64_t multi = 1;

    if(m < n)
    {
        quot = 0;
        *remain = m;
    }
    else if(m == n)
    {
        quot = 1;
        *remain = 0;
    }
    else
    {
        while(value >= n)
        {
            multi = 1;

            while((n*multi) <= (value>>1))
            {
                multi <<= 1;  //     2
            }
            quot += multi;
            value -= (n*multi);
        }
        *remain = value;
    }

    return quot;
}

uint64_t div3(uint64_t m, uint64_t n, uint64_t *remain)
{
    uint64_t quot = 0;
    uint64_t value = (m>>1);
    while(1)
    {
        if((n*(quot+value)) == m)
        {
            quot = quot+value;
            *remain = 0;
        }
        else if((n*(quot+value)) > m)
        {
            if(value == 1)
                break;
            else
                value >>= 1;
        }
        else
        {
            quot += value;

            if(value != 1)
                value >>= 1;
        }
    }
    *remain = m-(n*quot);
    return quot;

}

//     
int main(void)
{
    uint64_t m[4] = {0xffffffffffffffff, 261535774000000000, 60253400000, 4LL};
    uint64_t n[4] = {0xfffffffffffffffe, 6912470141400,      50698765,    2LL};
    int64_t quot = 0;        //   
    uint64_t remainder = 0;  //   
    int flag = 0;            //     
    int i = 0;


    for(i=1; i<4; i++)
    {
        //       
        if(m[i] < 0)
        {
            m[i] = -m[i];
            flag = 1;
        }
        if(n[i] < 0)
        {
            n[i] = -n[i];
            flag ^= 1;
        }

        quot = div1(m[i], n[i], &remainder);   //                

        if(flag)
            quot = -quot;

        printf("%d: %lld / %lld = %lld.....%lld
", i, m[i], n[i], quot, remainder);

        quot = 0;
        remainder = 0;
    }

    return 0;
}
【BZOJ 5152】【UOJ 347】【WC 2018】通路
負荷等化の基本アルゴリズム