0

wiki http://en.wikipedia.org/wiki/Karatsuba_algorithmの疑似コードに従っただけですが、 この実装の結果は非常に不安定です。時々動作しますが、100*100 のような場合です。失敗します。ここで見逃したものは何ですか？ご覧ください。

from math import *
f = lambda x: (int(x) & 1 and True) and 1
def fast_multiply( x = "100", y = "100"):
    print "input "+x+" | "+y
    int_buff = map( int, [x, y])
    if int_buff[0] < 10 or int_buff[1] < 10:
        #print "lol"
        return int_buff[0]*int_buff[1]

    degree = max( x.__len__(), y.__len__())

    higher_x, lower_x = x[ : int( ceil( len(x) / 2.0))], x[ len(x)/2 +f(len(x)):]
    higher_y, lower_y = y[ : int( ceil( len(y) / 2.0))], y[ len(y)/2 +f(len(y)):]
    #print lower_x+" & "+lower_y
    z0 = fast_multiply(lower_x, lower_y) #z0 = 0
    z1 = fast_multiply(str(int(lower_x)+int(higher_x)), str(int(lower_y)+int(higher_y)))
    z2 = fast_multiply(higher_x, higher_y)
    print "debug "+str(z0)+" "+str(z1)+" "+str(z2)
    return z2*(10**degree) + (z1-z2-z0)*(10**(degree/2))+z0




if __name__ == '__main__':
    print fast_multiply()

100*100 z2 が 100 になる場合は、これが正しいことに気付きました。これにより、 z2*(10**3)=100000 が得られますが、これは間違いです...

4

1 に答える 1

2

使用した疑似コードが間違っていました。問題はにありましたz2*(10**degree)。で計算するつもりだった2*m場所にベースを上げる必要がありました（両方とも である必要がありました）。mint( ceil(len(x) / 2.0))len(x)len(y)degree

私はそれをリファクタリングすることに抵抗できませんでした...少し。ウィキの定義からの名前を使用しました。任意の基数で実装するのは簡単ですが、簡単にするために 10 に固執しました。

def kmult(x, y):
    if min(x, y) < 10:
        return x * y

    m = half_ceil(degree(max(x, y)))

    x1, x0 = decompose(x, m)
    y1, y0 = decompose(y, m)

    z2 = kmult(x1, y1)
    z0 = kmult(x0, y0)
    z1 = kmult(x1 + x0, y1 + y0) - z2 - z0

    xy = z2 * 10**(2*m)  +  z1 * 10**m  +  z0
    return xy


def decompose(x, m):
    return x // 10 ** m, x % 10 ** m

def degree(x):
    return len(str(x))

def half_ceil(n):
    return n // 2 + (n & 1)

テスト:

print kmult(100, 100)

def test_kmult(r):
    for x, y in [(a, b) for b in range(r+1) for a in range(r+1)]:
        if kmult(x, y) != x * y:
            print('fail')
            break
    else:
        print('success')


test_kmult(100)

結果：

10000
success
于 2013-04-07T19:20:30.153 に答える