割り算が壊れたので自分で実装してみました。

あけましておめでとうございます。このブログは通常営業となりました。

概要

http://turi2.net/blog/724.html
http://d.hatena.ne.jp/nitoyon/20070629/four_operations_implementation_in_javascript
多倍長での除算がこの前うまく実装できなかったので、復讐のためやってみた。
ここではC言語で実装、比較や乗算などの使用も可能とする。ターゲットはunsigned int同士の除算で、unsigned long long intを使用した。

div1

順番に引いてくやつ。速度は結果に比例するので、普通は遅いが、このベンチマークでは結果の期待値が1なので速くでている。

div2

筆算。ビット数に比例する。このサイズなら速い。

div3

ニュートン法みたいな反復。ビット数の対数に比例する。このサイズだとdiv2より遅いが、多倍長なら話は別…のはず。

code

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>

/* slow method */
/* O(a/b) */
unsigned int div1(unsigned int a, unsigned int b)
{
    unsigned int ret = 0;
    while(a>=b) {
        ret++;
        a -= b;
    }
    return ret;
}

/* middle-speed method */
/* O(n). n is 32. */
unsigned int div2(unsigned int a, unsigned int b)
{
    unsigned int ret = 0;
    int i;
    for(i = __builtin_clz(b); i >= 0; i--) {
        if(a>=(b<<i)) {
            a -= b<<i;
            ret += 1<<i;
        }
    }
    return ret;
}

/* high-speed method */
/* O(log n). n is 32. */
unsigned int div3(unsigned int a, unsigned int b)
{
    int sb = 32-__builtin_clz(b);
    unsigned long long int ret = 1llu<<(32-sb);
    ret = (ret*((2llu<<32)-ret*b))>>32;
    ret = (ret*((2llu<<32)-ret*b))>>32;
    ret = (ret*((2llu<<32)-ret*b))>>32;
    ret = (ret*((2llu<<32)-ret*b))>>32;
    ret = (ret*((2llu<<32)-ret*b))>>32;
    ret = (ret*a)>>32;
    return (a-ret*b<b) ? ret : ret+1;
}
void checkdivs(unsigned int (*div)(unsigned int, unsigned int), int n) {
    struct timeval st, et;
    double ut;
    int i;
    int ccnt = 0;
    srand(10000);
    gettimeofday(&st, NULL);
    for(i = 0; i < n; i++) {
        unsigned int a = rand();
        unsigned int b = rand();
        if(b==0)b=1;
        if(a/b==div(a,b))ccnt++;
        else printf("%u / %u == %u != %u \n",a,b, a/b, div(a,b));
    }
    gettimeofday(&et, NULL);
    printf("%d/%d are correct.\n", ccnt, i);
    ut = (et.tv_usec-st.tv_usec)/1000.0 + (et.tv_sec-st.tv_sec)*1000.0;
    printf("time is %10.3f seconds.\n", ut/1000);
}

int main(int argc, char **argv)
{
    printf("div1: \n");checkdivs(div1, 10000000);
    printf("div2: \n");checkdivs(div2, 10000000);
    printf("div3: \n");checkdivs(div3, 10000000);
    return 0;
}

exec

div1:
10000000/10000000 are correct.
time is      1.201 seconds.
div2:
10000000/10000000 are correct.
time is      1.346 seconds.
div3:
10000000/10000000 are correct.
time is      1.944 seconds.