メタプログラミングによるべき和計算（2）

C++ Python

前回のべき和を求めるプログラムのアセンブリを見てみると、非常に複雑で時間がかかっていそうである。これなら、別のプログラムでべき和を求める公式を得たほうがいいのではないだろうか。それをPythonで求めてみた。



from itertools import *

from fractions import Fraction

import copy
def GCD(x, y):

    if x > y:

        return GCD(y, x)

    mod = y % x

    if mod == 0:

        return x

    else:

        return GCD(mod, x)
def LCM(x, y):

    return x * y / GCD(x, y)
def comb(n, m):

    if m == 0:

        return 1

    else:

        return comb(n, m - 1) * (n - m + 1) / m
def psub(x, y):

    z = copy.copy(x)

    for i in range(len(y)):

        z[i] -= y[i]

    return z
def pmul(x, y):

    z = [ 0 ] * (len(x) + len(y) - 1)

    for i in range(len(x)):

        for j in range(len(y)):

            z[i+j] += x[i] * y[j]

    return z
def pdiv(x, y):

    return map(lambda a: a * Fraction(1, y), x)
def ppow(x, n):

    return reduce(pmul, repeat(x, n), [ 1 ])
def S(n):

    p = ppow([1, 1], n + 1)

    p[0] -= 1

    for k in range(0, n):

        p = psub(p, pmul([ comb(n + 1, k) ], S(k)))

    return pdiv(p, n + 1)
def encode(a):

    N = reduce(LCM, map(lambda x: x.denominator, a), 1)

    b = map(lambda x: x * N, a)

    

    s = ""

    prev = 0

    for e in reversed(b):

        c = str(e)

        if s == "":

            s = c

        else:

            if prev == 0 or s.find('*') == -1:

                s += " * m"

            else:

                s = "(" + s + ") * m"

            if e < 0:

                s += " - " + str(-e)

            elif e > 0:

                s += " + " + c

        prev = e

    s += " / " + str(N)

    return s
a = S(5)

print encode(a)

前回と使う式は同じだが、それを使って求める公式の係数を計算している。最後にC++のコードに直している。出力は、



(((2 * m + 6) * m + 5) * m * m - 1) * m * m / 12

これをC++のコードに組み込んで、実行した。



#include 

#include 
#pragma comment( lib, "winmm.lib" )
using namespace std;
template

int S(int m) {

    return 1;

}
template<>

int S<5>(int m) {

    return (((2 * m + 6) * m + 5) * m * m - 1) * m * m / 12;

}
int main() {

    const int   N = 10000000;

    int sum = 0;

    int t = timeGetTime();

    for(int i = 0; i < N; i++) {

        int n = (rand() % 8 - 4);

        sum += S<5>(n);

    }

    cout << sum << endl;

    cout << (timeGetTime() - t) * 1e-3 << "s" << endl;

}

Sを前回のものに変えたものと比較すると、



前回：

774999921

2.087s
今回：

774999921

0.836s

それなりに速くなった。
ついでにSnの公式を20まで求めた。



S0 = m

S1 = (m2+m)/2

S2 = (2m3+3m2+m)/6

S3 = (m4+2m3+m2)/4

S4 = (6m5+15m4+10m3-m)/30

S5 = (2m6+6m5+5m4-m2)/12

S6 = (6m7+21m6+21m5-7m3+m)/42

S7 = (3m8+12m7+14m6-7m4+2m2)/24

S8 = (10m9+45m8+60m7-42m5+20m3-3m)/90

S9 = (2m10+10m9+15m8-14m6+10m4-3m2)/20

S10 = (6m11+33m10+55m9-66m7+66m5-33m3+5m)/66

S11 = (2m12+12m11+22m10-33m8+44m6-33m4+10m2)/24

S12 = (210m13+1365m12+2730m11-5005m9+8580m7-9009m5+4550m3-691m)/2730

S13 = (30m14+210m13+455m12-1001m10+2145m8-3003m6+2275m4-691m2)/420

S14 = (6m15+45m14+105m13-273m11+715m9-1287m7+1365m5-691m3+105m)/90

S15 = (3m16+24m15+60m14-182m12+572m10-1287m8+1820m6-1382m4+420m2)/48

S16 = (30m17+255m16+680m15-2380m13+8840m11-24310m9+44200m7-46988m5+23800m3-3617m)/510

S17 = (10m18+90m17+255m16-1020m14+4420m12-14586m10+33150m8-46988m6+35700m4-10851m2)/180

S18 = (210m19+1995m18+5985m17-27132m15+135660m13-529074m11+1469650m9-2678316m7+2848860m5-1443183m3+219335m)/3990

S19 = (42m20+420m19+1330m18-6783m16+38760m14-176358m12+587860m10-1339158m8+1899240m6-1443183m4+438670m2)/840

S20 = (330m21+3465m20+11550m19-65835m17+426360m15-2238390m13+8817900m11-24551230m9+44767800m7-47625039m5+24126850m3-3666831m)/6930