LoginSignup
38
25

More than 3 years have passed since last update.

@Red_Black_GPGPU

高速剰余変換で多倍長整数の乗算をGPUで実装して円周率3億桁計算した話

Last updated at Posted at 2019-10-23

0 はじめに

この記事では高速剰余変換と実装について解説していきます。この記事を書く上で参考にしたのはこのページです。
http://ushiro.jp/method/fmt.htm
http://www.cs.t-kougei.ac.jp/nsim/method/fmtbase.htm
基本的にはこのページで私が理解したとこ(+α)を解説します。多倍長整数の世界に興味のある方に少しでも面白さが伝わればと思い書きました。とか言ってますが私も全然素人です。ここではこの資料にのっとり高速剰余変換をFMT(Fast Modulo Transformation)と呼ぶこととします。 ←実は私がやってるのはNTT(数論変換)というようなものだったようで、詳しくは(FFTとNTTとFMTの違い)を見るとわかる通り剰余下でFFTやること＝FMTではなかったようで私はずっと勘違いしていました。以下FMTはNTTと適宜読み替えてください。

1 なぜ高速剰余変換(FMT)が必要か

モチベーションの根底に「多数桁の円周率をプログラムで求めたい」というのがあります。円周率自体は小数ですが、実際の計算はほとんど多倍長整数を使って計算します。ここで「多倍長整数の乗算をいかに高速に行うか」がプログラム全体の速度に大きく影響します。普通、N桁の乗算は愚直にやると$O(N^{2})$の計算量でとても遅いですが、FMTを使えばおおむね$O(NlogN)$のオーダーまで速くなります。
なおただ単に円周率を計算するならGMP等のライブラリやy-cruncherがあるのですが、最初からこれを使ってしまっては面白みがないので今は考えないこととします。

2 FMTの簡易実装例

なにとはともあれこれで10進4桁の整数の乗算をしてみましょう。
Python3の簡易版FMT関数と逆FMT関数です。

fmt_smaple0.py
import numpy as np

# MODP剰余下でaのb乗を返す
# 追記1:pow(a,b,MODP)を使ったほうが速い!
def ModExp(a, b):
    ans = np.uint64(1)
    aa = np.uint64(a)
    while (b != 0):
        if (b % 2 == 1):
            ans = ans * aa % np.uint64(MODP)
        aa = aa * aa % np.uint64(MODP)
        b //= 2
    return np.uint32(ans)

# bit逆順 uint型を想定
def bitRev(a):
    a = (a & 0x55555555) << 1 | (a & 0xAAAAAAAA) >> 1
    a = (a & 0x33333333) << 2 | (a & 0xCCCCCCCC) >> 2
    a = (a & 0x0F0F0F0F) << 4 | (a & 0xF0F0F0F0) >> 4
    a = (a & 0x00FF00FF) << 8 | (a & 0xFF00FF00) >> 8
    a = (a & 0x0000FFFF) << 16 | (a & 0xFFFF0000) >> 16
    return a >> (32 - digit_level)

def FMT(inputA):
    # まず配列をビット逆順
    arrayA = inputA[bitRev(np.uint32(range(digitN)))]
    # FMTループ
    for i in range(digit_level):
        loopCnt_Pow2 = 1 << i
        for idx in range(digitN // 2):
            t2 = idx % loopCnt_Pow2
            t0 = idx * 2 - t2
            t1 = t0 + loopCnt_Pow2
            arrayAt0 = arrayA[t0]
            arrayAt1 = arrayA[t1]
            w0 = ModExp(MODP_Wn, t2 * (1 << (digit_level - 1 - i)))
            w1 = np.uint32(np.uint64(arrayAt1) * np.uint64(w0) % np.uint64(MODP))
            arrayA[t1] = (arrayAt0 + MODP - w1) % MODP
            arrayA[t0] = (arrayAt0 + w1) % MODP
    return arrayA[:]

def iFMT(inputA):
    # まず配列をビット逆順
    arrayA = inputA[bitRev(np.uint32(range(digitN)))]
    # FMTループ
    for i in range(digit_level):
        loopCnt_Pow2 = 1 << i
        for idx in range(digitN // 2):
            t2 = idx % loopCnt_Pow2
            t0 = idx * 2 - t2
            t1 = t0 + loopCnt_Pow2
            arrayAt0 = arrayA[t0]
            arrayAt1 = arrayA[t1]
            w0 = ModExp(MODP_Wn, digitN - t2 * (1 << (digit_level - 1 - i)))
            w1 = np.uint32(np.uint64(arrayAt1) * np.uint64(w0) % np.uint64(MODP))
            arrayA[t1] = (arrayAt0 + MODP - w1) % MODP
            arrayA[t0] = (arrayAt0 + w1) % MODP
    return arrayA[:]

def DivN(inputA):
    digitN_r = ModExp(digitN, MODP - 2)  # digitNのモジュラ逆数を計算
    ret = np.uint64(inputA[:]) * np.uint64(digitN_r) % np.uint64(MODP)
    return np.uint32(ret)


if __name__=="__main__":
    # 可変パラメーター
    digit_level = 3
    digitN = 1 << digit_level
    MODP = np.uint32(257)
    MODP_Wn = np.uint32(4) # 1の原始N乗根でなければならない。4の8乗=65536。65536%257=1

    # C=A*Bをしたい
    A=np.uint32([1,4,1,4,0,0,0,0])
    B=np.uint32([2,1,3,5,0,0,0,0])

    A = FMT(A)  # 変換
    B = FMT(B)  # 変換
    C = np.uint64(A[:]) * np.uint64(B[:]) % np.uint64(MODP)  # 各要素どうし乗算
    C = np.uint32(C)
    C = iFMT(C)  # 逆変換
    C = DivN(C)  # Nで割る

    print("畳み込みの結果")
    for i in range(digitN):
        print("C[",i,"] =",C[i])

A,B各要素に10進1桁が入っています。ここでは
A=np.uint32([1,4,1,4,0,0,0,0])
B=np.uint32([2,1,3,5,0,0,0,0])
としているので4141*5312の計算をしようとしています。答えは21996992となるはずです。

実行結果
畳み込みの結果
C[ 0 ] = 2
C[ 1 ] = 9
C[ 2 ] = 9
C[ 3 ] = 26
C[ 4 ] = 27
C[ 5 ] = 17
C[ 6 ] = 20
C[ 7 ] = 0

C[0]が1の位、C[1]が10の位・・・で畳み込みの結果が出力されています。
まだ繰り上がりの計算をしていませんので、計算すると

実行結果
C[ 0 ] = 2
C[ 1 ] = 9
C[ 2 ] = 9
C[ 3 ] = 6
C[ 4 ] = 9
C[ 5 ] = 9
C[ 6 ] = 1
C[ 7 ] = 2

2・1・9・9・6・9・9・2
あってます！

次は1234*5678をやってみます。答えは7006652になるはずです。

実行結果
畳み込みの結果
C[ 0 ] = 32
C[ 1 ] = 52
C[ 2 ] = 61
C[ 3 ] = 60
C[ 4 ] = 34
C[ 5 ] = 16
C[ 6 ] = 5
C[ 7 ] = 0
繰り上がり処理の結果
C[ 0 ] = 2
C[ 1 ] = 5
C[ 2 ] = 6
C[ 3 ] = 6
C[ 4 ] = 0
C[ 5 ] = 0
C[ 6 ] = 7
C[ 7 ] = 0

7・0・0・6・6・5・2
あってます！
面白いですね。

次は9999*9999をやってみます。

実行結果
畳み込みの結果
C[ 0 ] = 81
C[ 1 ] = 162
C[ 2 ] = 243
C[ 3 ] = 67
C[ 4 ] = 243
C[ 5 ] = 162
C[ 6 ] = 81
C[ 7 ] = 0
繰り上がり処理の結果
C[ 0 ] = 1
C[ 1 ] = 0
C[ 2 ] = 0
C[ 3 ] = 3
C[ 4 ] = 2
C[ 5 ] = 7
C[ 6 ] = 9
C[ 7 ] = 9

答えは99980001になるはずですが、残念ながら違った結果がでています。
この簡易版FMTでは法$P=257$としており、すべての計算は257の剰余下で行われています。本当は畳み込み結果C[3]が324のところを67にされているのが原因です。剰余版オーバーフローみたいなもんです。

ここで少しFMT関数が何をしているかを説明します。
FMT関数はaを入力としてこのようなfを出力しています。

$f_k=a_0ω^{0}+a_1ω^{k}+a_2ω^{2k}....+a_{n-1}ω^{(n-1)k}$

(k=0,1,2・・・・n-1)

$ω$は適当な整数とし、$ω^{n}≡1$ $(mod$ $P)$となるPの剰余下で計算をします。$ω$と$P$は互いに素な必要があります。これにより$ω^1,ω^2,ω^3・・$は$ω^n$まで同じ数が出現せず、$ω^n$でちょうど1周回ってくる性質のためFFTと同じくバラフライ演算ができ、FFTと同様高速に計算できるのです。逆変換は$ω$を$-ω$に置き換えるだけです(ただし最後全体をNで割るのはFFTと同じ)。
FFTでは小数や複素数で計算を進めていきますが、これを全て整数でやることで小数特有の誤差を気にする必要がなくなります。ただ先ほど述べたように全ての計算を整数剰余下で行うので音声信号の周波数を求めたりする用途等には使われず、多倍長整数の乗算や競プロの問題に利用される程度で、FFTと比べたら用途は狭いようです。
多倍長整数の乗算を行うためには上記のコードのように

A,BをFMTで変換
↓
変換後のA,Bの各要素同士を乗算
↓
FMT逆変換

をすることで、畳み込み結果を得ることができます。

畳み込み結果もmod Pされているので、畳み込み結果がP以上になる可能性がある場合、計算が正しく行えません。なのでP自体を大きくするか、複数のPを使って剰余下の結果を得た後、中国余剰定理で復元するしかありません。基本的には後者のアプローチをとることが多いです。

参考
【FFT,FMTで$O(NlogN$)で畳み込みできる原理】
https://qiita.com/peria/items/fbdd52768b4659823d88

2-1 Pを大きくとるアプローチ

とりあえず257は小さすぎるので$P=65537$でやってみます。ソースコードのMODPとMODP_Wnを以下に書き換えます。
MODP = 65537
MODP_Wn = 16
これでやると

畳み込みの結果
C[ 0 ] = 81
C[ 1 ] = 162
C[ 2 ] = 243
C[ 3 ] = 324
C[ 4 ] = 243
C[ 5 ] = 162
C[ 6 ] = 81
C[ 7 ] = 0
繰り上がり処理の結果
C[ 0 ] = 1
C[ 1 ] = 0
C[ 2 ] = 0
C[ 3 ] = 0
C[ 4 ] = 8
C[ 5 ] = 9
C[ 6 ] = 9
C[ 7 ] = 9

99980001
あってます。
しかし、ただPを大きくするのでは限界があります。もっと大きい桁数で計算することを考えるとPの必要最低桁数は2進64桁以上になりそうで、いよいよ剰余計算自体がlong型でも扱えなくなってしまいます。

2-2 中国余剰定理を使ったアプローチ

Pを257と17でやってみましょう

fmt_smaple1.py
import numpy as np

# MODP剰余下でaのb乗を返す
# 追記1:np.uint32(pow(np.uint32(a),np.uint32(b),np.uint32(MODP)))を使ったほうが速い!
def ModExp(a, b):
    ans = np.uint64(1)
    aa = np.uint64(a)
    while (b != 0):
        if (b % 2 == 1):
            ans = ans * aa % np.uint64(MODP)
        aa = aa * aa % np.uint64(MODP)
        b //= 2
    return np.uint32(ans)

# aの逆元を返す mod m
def modinv(a,m):
    b = m
    u = 1
    v = 0
    while (b) :
        t = a // b
        a -= t * b
        tmp=a
        a=b
        b=tmp
        u -= t * v
        tmp = u
        u = v
        v = tmp
    u %= m
    if u < 0:
        u += m
    return u

# bit逆順 uint型を想定
def bitRev(a):
    a = (a & 0x55555555) << 1 | (a & 0xAAAAAAAA) >> 1
    a = (a & 0x33333333) << 2 | (a & 0xCCCCCCCC) >> 2
    a = (a & 0x0F0F0F0F) << 4 | (a & 0xF0F0F0F0) >> 4
    a = (a & 0x00FF00FF) << 8 | (a & 0xFF00FF00) >> 8
    a = (a & 0x0000FFFF) << 16 | (a & 0xFFFF0000) >> 16
    return a >> (32 - digit_level)

def FMT(inputA):
    # まず配列をビット逆順
    arrayA = inputA[bitRev(np.uint32(range(digitN)))]
    # FMTループ
    for i in range(digit_level):
        loopCnt_Pow2 = 1 << i
        for idx in range(digitN // 2):
            t2 = idx % loopCnt_Pow2
            t0 = idx * 2 - t2
            t1 = t0 + loopCnt_Pow2
            arrayAt0 = arrayA[t0]
            arrayAt1 = arrayA[t1]
            w0 = ModExp(MODP_Wn, t2 * (1 << (digit_level - 1 - i)))
            w1 = np.uint32(np.uint64(arrayAt1) * np.uint64(w0) % np.uint64(MODP))
            arrayA[t1] = (arrayAt0 + MODP - w1) % MODP
            arrayA[t0] = (arrayAt0 + w1) % MODP
    return arrayA[:]

def iFMT(inputA):
    # まず配列をビット逆順
    arrayA = inputA[bitRev(np.uint32(range(digitN)))]
    # FMTループ
    for i in range(digit_level):
        loopCnt_Pow2 = 1 << i
        for idx in range(digitN // 2):
            t2 = idx % loopCnt_Pow2
            t0 = idx * 2 - t2
            t1 = t0 + loopCnt_Pow2
            arrayAt0 = arrayA[t0]
            arrayAt1 = arrayA[t1]
            w0 = ModExp(MODP_Wn, digitN - t2 * (1 << (digit_level - 1 - i)))
            w1 = np.uint32(np.uint64(arrayAt1) * np.uint64(w0) % np.uint64(MODP))
            arrayA[t1] = (arrayAt0 + MODP - w1) % MODP
            arrayA[t0] = (arrayAt0 + w1) % MODP
    return arrayA[:]

def DivN(inputA):
    digitN_r = modinv(digitN, int(MODP))  # digitNの逆元を計算
    ret = np.uint64(inputA[:]) * np.uint64(digitN_r) % np.uint64(MODP)
    return np.uint32(ret)


if __name__=="__main__":
    digit_level = 3
    digitN = 1 << digit_level

    # C=A*Bをしたい
    A_=np.uint32([9,9,9,9,0,0,0,0])
    B_=np.uint32([9,9,9,9,0,0,0,0])

    #まずp=257でやる
    MODP = np.uint32(257)
    MODP_Wn = np.uint32(4) # 1の原始N乗根でなければならない。4の8乗=65536。65536%257=1

    A = FMT(A_)  # 変換
    B = FMT(B_)  # 変換
    C0 = np.uint64(A[:]) * np.uint64(B[:]) % np.uint64(MODP)  # 各要素どうし乗算
    C0 = np.uint32(C0)
    C0 = iFMT(C0)  # 逆変換
    C0 = DivN(C0)  # Nで割る

    #p=17でやる
    MODP = np.uint32(17)
    MODP_Wn = np.uint32(2) # 1の原始N乗根でなければならない。2の8乗=256。256%17=1

    A = FMT(A_)  # 変換
    B = FMT(B_)  # 変換
    C1 = np.uint64(A[:]) * np.uint64(B[:]) % np.uint64(MODP)  # 各要素どうし乗算
    C1 = iFMT(C1)  # 逆変換
    C1 = DivN(C1)  # Nで割る

    print("畳み込みの結果 mod","257","\t\t畳み込みの結果 mod","17","\t\tGarnerで復元")
    for i in range(digitN):
        print("C[",i,"] =",C0[i],
              " \t\t\t\t","C[",i,"] =",C1[i],
              "\t\t\t\t",int(C0[i])+((int(C1[i])-int(C0[i]))*modinv(257,17))%17*257
              )

実行結果↓

畳み込みの結果 mod 257       畳み込みの結果 mod 17        Garnerで復元
C[ 0 ] = 81                  C[ 0 ] = 13                 81
C[ 1 ] = 162                 C[ 1 ] = 9                  162
C[ 2 ] = 243                 C[ 2 ] = 5                  243
C[ 3 ] = 67                  C[ 3 ] = 1                  324
C[ 4 ] = 243                 C[ 4 ] = 5                  243
C[ 5 ] = 162                 C[ 5 ] = 9                  162
C[ 6 ] = 81                  C[ 6 ] = 13                 81
C[ 7 ] = 0                   C[ 7 ] = 0                  0

ここではGarnerのアルゴリズムというのを使って復元しています。ちゃんとC[3]が324になってますね！復元したものを繰り上がり処理すれば
99980001
が得られます。

Garnerアルゴリズムについて超ざっくり説明すると、復元したいxがあって
x%a が ar
x%b が br
のとき

ar+(br-ar)*modinv(a,b)%b*a

でもとの値が復元できるというアルゴリズムです。modinvは「aをかけるとbの剰余下で1になる整数」を返す関数です。
Garnerのアルゴリズムはわかりやすいサイトがいくつもあるのでそちらをご参照ください。
https://www.creativ.xyz/ect-gcd-crt-garner-927/
https://qiita.com/drken/items/ae02240cd1f8edfc86fd

他にも拡張ユークリッド互除法で値を復元することもでき、それを使ってもいいと思います。

このようにPの値自体を大きくせずとも、使うPの種類を増やしていけば大きな数を復元できることがわかります。

3 大きめの桁に挑戦

3-1 1要素に32bit敷き詰める

今までわかりやすさ優先で10進法前提でやっていましたが4294967296進法(=2の32乗)で多倍長乗算をやりたいと思います。

    A=np.uint32([4294967295,4294967295,4294967295,4294967295,0,0,0,0])
    B=np.uint32([4294967295,4294967295,4294967295,4294967295,0,0,0,0])

入力側を32bit敷き詰めました。右半分はまだ0ですがこれもあとで使えるようになります。
このまま動かすとGarnerで復元できる257*17を大幅に超えるので法Pも大きくします。
MODPを469762049 , 1811939329 , 2013265921
MODP_Wnを443138433 , 1452317833 , 1801542727
とします。
あと「FMT→同じ要素同士乗算→iFMT→Nで割る」の一連の操作をComvolution関数にまとめて全体をクラス化します。

fmt_smaple2.py
import numpy as np

# aの逆元を返す mod m
def modinv(a,m):
    b = m
    u = 1
    v = 0
    while (b) :
        t = a // b
        a -= t * b
        tmp=a
        a=b
        b=tmp
        u -= t * v
        tmp = u
        u = v
        v = tmp
    u %= m
    if u < 0:
        u += m
    return u

# bit逆順 uint型を想定
def bitRev(a):
    a = (a & 0x55555555) << 1 | (a & 0xAAAAAAAA) >> 1
    a = (a & 0x33333333) << 2 | (a & 0xCCCCCCCC) >> 2
    a = (a & 0x0F0F0F0F) << 4 | (a & 0xF0F0F0F0) >> 4
    a = (a & 0x00FF00FF) << 8 | (a & 0xFF00FF00) >> 8
    a = (a & 0x0000FFFF) << 16 | (a & 0xFFFF0000) >> 16
    return a >> (32 - digit_level)

#garnerのアルゴリズム、3変数ver
#64bit整数におさまらない可能性あり、pythonのint型を使用
def Garner(a_,b_,c_,ar_,br_,cr_):
    a = int(a_)
    b = int(b_)
    c = int(c_)
    ar = int(ar_)
    br = int(br_)
    cr = int(cr_)
    x = ar + (br - ar) * modinv(a,b) % b * a
    x = x + (cr - x % c) * modinv(a,c) % c * modinv(b,c) % c * a * b
    #x %= a*b*c #a<b<cのときこの処理は必要ない
    return x

class FMTClass():
    def __init__(self,MODP,MODP_Wn):
        self.MODP=np.uint32(MODP)
        self.MODP_Wn=np.uint32(MODP_Wn)

    # MODP剰余下でaのb乗を返す
    # 追記1:np.uint32(pow(np.uint32(a),np.uint32(b),np.uint32(self.MODP)))を使ったほうが速い!
    def ModExp(self,a, b):
        ans = np.uint64(1)
        aa = np.uint64(a)
        while (b != 0):
            if (b % 2 == 1):
                ans = ans * aa % np.uint64(self.MODP)
            aa = aa * aa % np.uint64(self.MODP)
            b //= 2
        return np.uint32(ans)

    def FMT(self,inputA):
        # まず配列をビット逆順
        arrayA = inputA[bitRev(np.uint32(range(digitN)))]
        # FMTループ
        for i in range(digit_level):
            loopCnt_Pow2 = 1 << i
            for idx in range(digitN // 2):
                t2 = idx % loopCnt_Pow2
                t0 = idx * 2 - t2
                t1 = t0 + loopCnt_Pow2
                arrayAt0 = arrayA[t0]
                arrayAt1 = arrayA[t1]
                w0 = self.ModExp(self.MODP_Wn, t2 * (1 << (digit_level - 1 - i)))
                w1 = np.uint32(np.uint64(arrayAt1) * np.uint64(w0) % np.uint64(self.MODP))
                arrayA[t1] = (arrayAt0 + self.MODP - w1) % self.MODP
                arrayA[t0] = (arrayAt0 + w1) % self.MODP
        return arrayA[:]

    def iFMT(self,inputA):
        # まず配列をビット逆順
        arrayA = inputA[bitRev(np.uint32(range(digitN)))]
        # FMTループ
        for i in range(digit_level):
            loopCnt_Pow2 = 1 << i
            for idx in range(digitN // 2):
                t2 = idx % loopCnt_Pow2
                t0 = idx * 2 - t2
                t1 = t0 + loopCnt_Pow2
                arrayAt0 = arrayA[t0]
                arrayAt1 = arrayA[t1]
                w0 = self.ModExp(self.MODP_Wn, digitN - t2 * (1 << (digit_level - 1 - i)))
                w1 = np.uint32(np.uint64(arrayAt1) * np.uint64(w0) % np.uint64(self.MODP))
                arrayA[t1] = (arrayAt0 + self.MODP - w1) % self.MODP
                arrayA[t0] = (arrayAt0 + w1) % self.MODP
        return arrayA[:]

    def DivN(self,inputA):
        digitN_r = modinv(digitN, int(self.MODP))  # digitNの逆元を計算
        ret = np.uint64(inputA[:]) * np.uint64(digitN_r) % np.uint64(self.MODP)
        return np.uint32(ret)

    #A*Bの畳み込み結果を返す
    def Convolution(self,A_,B_):
        A = self.FMT(A_%self.MODP)
        B = self.FMT(B_%self.MODP)
        C = np.uint64(A[:]) * np.uint64(B[:]) % np.uint64(self.MODP)  # 各要素どうし乗算
        C = np.uint32(C) #結果は32bit範囲内におさまってるので
        C = self.iFMT(C)  # 逆変換
        C = self.DivN(C)  # Nで割る
        return C



if __name__=="__main__":
    digit_level = 3
    digitN = 1 << digit_level

    FMT_P0 = FMTClass(MODP=469762049 , MODP_Wn=443138433)
    FMT_P1 = FMTClass(MODP=1811939329, MODP_Wn=1452317833)
    FMT_P2 = FMTClass(MODP=2013265921, MODP_Wn=1801542727)

    # C=A*Bをしたい
    A_=np.uint32([4294967295,4294967295,4294967295,4294967295,0,0,0,0])
    B_=np.uint32([4294967295,4294967295,4294967295,4294967295,0,0,0,0])

    #p=469762049でやる
    C0 = FMT_P0.Convolution(A_,B_)

    #p=1811939329でやる
    C1 = FMT_P1.Convolution(A_, B_)

    #p=2013265921でやる
    C2 = FMT_P2.Convolution(A_, B_)

    print("畳み込みの結果 mod","469762049",
          "\t\t畳み込みの結果 mod","1811939329",
          "\t\t畳み込みの結果 mod", "2013265921",
          "\t\tGarnerで復元")

    for i in range(digitN):
        print("C[",i,"] =",C0[i],
              "  \t\t\t\t", "C[", i, "] =", C1[i],
              " \t\t\t\t", "C[", i, "] =", C2[i],
              "\t\t\t\t",Garner(FMT_P0.MODP,FMT_P1.MODP,FMT_P2.MODP,C0[i],C1[i],C2[i])
              )

結果は

畳み込みの結果 mod 469762049         畳み込みの結果 mod 1811939329        畳み込みの結果 mod 2013265921        Garnerで復元
C[ 0 ] = 325957443                   C[ 0 ] = 1429170264                 C[ 0 ] = 635297256                  18446744065119617025
C[ 1 ] = 182152837                   C[ 1 ] = 1046401199                 C[ 1 ] = 1270594512                 36893488130239234050
C[ 2 ] = 38348231                    C[ 2 ] = 663632134                  C[ 2 ] = 1905891768                 55340232195358851075
C[ 3 ] = 364305674                   C[ 3 ] = 280863069                  C[ 3 ] = 527923103                  73786976260478468100
C[ 4 ] = 38348231                    C[ 4 ] = 663632134                  C[ 4 ] = 1905891768                 55340232195358851075
C[ 5 ] = 182152837                   C[ 5 ] = 1046401199                 C[ 5 ] = 1270594512                 36893488130239234050
C[ 6 ] = 325957443                   C[ 6 ] = 1429170264                 C[ 6 ] = 635297256                  18446744065119617025
C[ 7 ] = 0                   C[ 7 ] = 0                  C[ 7 ] = 0                  0

検算してみましょう。

#愚直に計算
import numpy as np
A_ = np.uint32([4294967295, 4294967295, 4294967295, 4294967295, 0, 0, 0, 0])
B_ = np.uint32([4294967295, 4294967295, 4294967295, 4294967295, 0, 0, 0, 0])

A=0
B=0
r=1
for i in range(8):
    A += int(A_[i]) * r
    B += int(B_[i]) * r
    r*=2**32
print(A*B)

#FMTの答え
C_=[18446744065119617025,36893488130239234050,55340232195358851075,73786976260478468100,55340232195358851075,36893488130239234050,18446744065119617025,0]
C=0
r=1
for i in range(8):
    C+=C_[i]*r
    r*=2**32
print(C)

結果は

115792089237316195423570985008687907852589419931798687112530834793049593217025
115792089237316195423570985008687907852589419931798687112530834793049593217025

あってました！嬉しいですね！
これで4294967296進数で多倍長整数の乗算ができるようになりました。np.uint32型に32bit敷き詰めることができて気持ちいいです。また、Garnerで復元するところ以外は、計算は全て64bit整数の範囲内におさまっています。Garnerで復元するところはどうしても64bit整数をオーバーするので「Python3のint型」を使っています(Python3のint型は実質多倍長整数)。ここはGPU(CUDA)に移植する際には少し考えないといけなさそうですね。

3-2 A,Bの要素数を増やす

今までわかりやすさ優先でAとBを8要素でやっていましたが、もっと増やしてみたいと思います。要素数が$2^N$ならなんでも良い・・・と言いたいところですが、前述のよう$ω^n≡1$ $(mod$ $P)$という条件があります。

先ほどMODPを$469762049$ , $1811939329$ , $2013265921$
MODP_Wnを$443138433$ , $1452317833$ , $1801542727$
としていましたが、それぞれ

$443138433^8≡1$ $(mod$ $469762049)$
$1452317833^8≡1$ $(mod$ $1811939329)$
$1801542727^8≡1$ $(mod$ $2013265921)$

が成りたっていました。
$443138433$は$8$乗すると$469762049$剰余下で$1$になるので$443138433$は$1$の原始$8$乗根と言います。

要素数を$2^N$にしたいのであれば、1の原始$2^N$乗根とPを用意しなければいけませんが、ここでその例をいくつかあげてみます。

$52278$は$2013265921$を法としたとき、$1$の原始$2^{27}$乗根
$60733$は$469762049$を法としたとき、$1$の原始$2^{26}$乗根
$59189$は$1811939329$を法としたとき、$1$の原始$2^{26}$乗根
です。
ちょうどいいのでこれを使います。これでいっきに要素数を$2^{26}$まで増やして、1億桁近くの多倍長乗算ができます！

と、思ってやってみたらさすがに重くて結果が出るまで2日かかりました。やっぱりPythonは遅いですね()。数十秒で終わるように要素数を$2^{14}$に抑えたバージョンを貼ります。

fmt_smaple3.py
import numpy as np

# aの逆元を返す mod m
def modinv(a,m):
    b = m
    u = 1
    v = 0
    while (b) :
        t = a // b
        a -= t * b
        tmp=a
        a=b
        b=tmp
        u -= t * v
        tmp = u
        u = v
        v = tmp
    u %= m
    if u < 0:
        u += m
    return u

# bit逆順 uint型を想定
def bitRev(a):
    a = (a & 0x55555555) << 1 | (a & 0xAAAAAAAA) >> 1
    a = (a & 0x33333333) << 2 | (a & 0xCCCCCCCC) >> 2
    a = (a & 0x0F0F0F0F) << 4 | (a & 0xF0F0F0F0) >> 4
    a = (a & 0x00FF00FF) << 8 | (a & 0xFF00FF00) >> 8
    a = (a & 0x0000FFFF) << 16 | (a & 0xFFFF0000) >> 16
    return a >> (32 - digit_level)

#garnerのアルゴリズム、3変数ver
#64bit整数におさまらない可能性あり、pythonのint型を使用
def Garner(a_,b_,c_,ar_,br_,cr_):
    a = int(a_)
    b = int(b_)
    c = int(c_)
    ar = int(ar_)
    br = int(br_)
    cr = int(cr_)
    x = ar + (br - ar) * modinv(a,b) % b * a
    x = x + (cr - x % c) * modinv(a,c) % c * modinv(b,c) % c * a * b
    #x %= a*b*c #a<b<cのときこの処理は必要ない
    return x

class FMTClass():
    def __init__(self,MODP,MODP_Wn):
        self.MODP=np.uint32(MODP)
        self.MODP_Wn=np.uint32(MODP_Wn)

    # MODP剰余下でaのb乗を返す
    # 追記1:np.uint32(pow(np.uint32(a),np.uint32(b),np.uint32(self.MODP)))を使ったほうが速い!
    def ModExp(self,a, b):
        ans = np.uint64(1)
        aa = np.uint64(a)
        while (b != 0):
            if (b % 2 == 1):
                ans = ans * aa % np.uint64(self.MODP)
            aa = aa * aa % np.uint64(self.MODP)
            b //= 2
        return np.uint32(ans)

    def FMT(self,inputA):
        # まず配列をビット逆順
        arrayA = inputA[bitRev(np.uint32(range(digitN)))]
        # FMTループ
        for i in range(digit_level):
            loopCnt_Pow2 = 1 << i
            for idx in range(digitN // 2):
                t2 = idx % loopCnt_Pow2
                t0 = idx * 2 - t2
                t1 = t0 + loopCnt_Pow2
                arrayAt0 = arrayA[t0]
                arrayAt1 = arrayA[t1]
                w0 = self.ModExp(self.MODP_Wn, t2 * (1 << (digit_level - 1 - i)))
                w1 = np.uint32(np.uint64(arrayAt1) * np.uint64(w0) % np.uint64(self.MODP))
                arrayA[t1] = (arrayAt0 + self.MODP - w1) % self.MODP
                arrayA[t0] = (arrayAt0 + w1) % self.MODP
        return arrayA[:]

    def iFMT(self,inputA):
        # まず配列をビット逆順
        arrayA = inputA[bitRev(np.uint32(range(digitN)))]
        # FMTループ
        for i in range(digit_level):
            loopCnt_Pow2 = 1 << i
            for idx in range(digitN // 2):
                t2 = idx % loopCnt_Pow2
                t0 = idx * 2 - t2
                t1 = t0 + loopCnt_Pow2
                arrayAt0 = arrayA[t0]
                arrayAt1 = arrayA[t1]
                w0 = self.ModExp(self.MODP_Wn, digitN - t2 * (1 << (digit_level - 1 - i)))
                w1 = np.uint32(np.uint64(arrayAt1) * np.uint64(w0) % np.uint64(self.MODP))
                arrayA[t1] = (arrayAt0 + self.MODP - w1) % self.MODP
                arrayA[t0] = (arrayAt0 + w1) % self.MODP
        return arrayA[:]

    def DivN(self,inputA):
        digitN_r = modinv(digitN, int(self.MODP))  # digitNの逆元を計算
        ret = np.uint64(inputA[:]) * np.uint64(digitN_r) % np.uint64(self.MODP)
        return np.uint32(ret)

    #A*Bの畳み込み結果を返す
    def Convolution(self,A_,B_):
        A = self.FMT(A_%self.MODP)
        B = self.FMT(B_%self.MODP)
        C = np.uint64(A[:]) * np.uint64(B[:]) % np.uint64(self.MODP)  # 各要素どうし乗算
        C = np.uint32(C) #結果は32bit範囲内におさまってるので
        C = self.iFMT(C)  # 逆変換
        C = self.DivN(C)  # Nで割る
        return C

#検算用。A,Bの配列から一つの多倍長整数を作る。pythonのint型は乗算にKaratsuba法を使っているのでちょっと速い
def CalcAB(A,B):
    Alist = A.tolist()
    Blist = B.tolist()
    r=2**32
    for i in range(digit_level):
        for j in range(2**(digit_level-1-i)):
            Alist[j] = Alist[j * 2] + Alist[j * 2 + 1] * r
            Blist[j] = Blist[j * 2] + Blist[j * 2 + 1] * r
        r = r * r
    return Alist[0] * Blist[0]

#検算用。Cから一つの多倍長整数を作る
def CalcC(C):
    Clist = C.tolist()
    r = 2 ** 32
    for i in range(digit_level):
        for j in range(2 ** (digit_level-1 - i)):
            Clist[j] = Clist[j * 2] + Clist[j * 2 + 1] * r
        r = r * r
    return Clist[0]

def CreateWn(a,m,n_level):
    ans=a
    for i in range(n_level-digit_level):
        ans=ans*ans%m
    return ans

#E0,E1,E2を使ってgarnerでEを復元して繰り上がり処理
def Carrying(C0,C1,C2):
    C = np.zeros(digitN).astype(np.uint32)
    upg = 0  # 繰り上がり分 python3のint型なのでオーバーフローの心配はない
    for i in range(digitN):
        g = Garner(FMT_P0.MODP, FMT_P1.MODP, FMT_P2.MODP, C0[i], C1[i], C2[i])
        g += upg
        C[i] = np.uint32(g % 4294967296)
        upg = g//4294967296
    return C


if __name__=="__main__":
    digit_level = 14
    digitN = 1 << digit_level

    FMT_P0 = FMTClass(MODP=469762049 , MODP_Wn=CreateWn(60733,469762049,26))
    FMT_P1 = FMTClass(MODP=1811939329, MODP_Wn=CreateWn(59189,1811939329,26))
    FMT_P2 = FMTClass(MODP=2013265921, MODP_Wn=CreateWn(52278,2013265921,27))

    # C=A*Bをしたい
    # A,Bの初期化。前半にランダムな値をいれる。あとは0
    A_ = np.zeros(digitN).astype(np.uint32)
    A_[0:digitN // 2] = np.uint32( np.random.randint(0, 2147483648, (digitN // 2)) * 2
                                   + np.random.randint(0, 2, (digitN // 2)))
    B_ = np.zeros(digitN).astype(np.uint32)
    B_[0:digitN // 2] = np.uint32(np.random.randint(0, 2147483648, (digitN // 2)) * 2
                                  + np.random.randint(0, 2, (digitN // 2)))

    #p=469762049でやる
    C0 = FMT_P0.Convolution(A_,B_)

    #p=1811939329でやる
    C1 = FMT_P1.Convolution(A_, B_)

    #p=2013265921でやる
    C2 = FMT_P2.Convolution(A_, B_)

    #一部結果表示
    print("畳み込みの結果 mod","469762049",
          "\t\t畳み込みの結果 mod","1811939329",
          "\t\t畳み込みの結果 mod", "2013265921",
          "\t\tGarnerで復元")
    for i in range(10):
        print("C[",i,"] =",C0[i],
              "  \t\t\t\t", "C[", i, "] =", C1[i],
              " \t\t\t\t", "C[", i, "] =", C2[i],
              "\t\t\t\t",Garner(FMT_P0.MODP,FMT_P1.MODP,FMT_P2.MODP,C0[i],C1[i],C2[i])
              )

    #繰り上がり
    C=Carrying(C0,C1,C2)

    #検算
    print("検算開始")
    if CalcAB(A_, B_)==CalcC(C):
        print("正解！")
    else:
        print("ちがう！")

なんだか検算のほうが速かった気がします・・・。(まぁとりあえず今は正しく動いていることが確認できただけでヨシとします。)

上記コードですが


FMT_P0 = FMTClass(MODP=469762049 , MODP_Wn=CreateWn(60733,469762049,26))

このインスタンスに初期値を設定するところで、CreateWn関数を使っています。ここで要素数$2^N$に応じた原始$2^N$乗根を生成している感じです。例えば$60733$は$1$の$2^{26}$乗根と分かっているので、$2^{25}$乗根は$60733*60733$ $(mod$ $469762049)$で簡単に計算できます。

また余談ですが、自前で$P$とそれに対応する原始$2^N$乗根を探しあてるには全探索しても良いのですが、フェルマーの小定理を知っておくと良いかもしれません。今回使っているPはそれぞれ
$2013265921=15×2^{27}+1$
$1811939329=27×2^{26}+1$
$469762049 = 7×2^{26}+1$
であり、それぞれ$2^{27}$乗根、$2^{26}$乗根、$2^{26}$乗根まで存在します。
私は$2^N$に適当な数をかけて1を足した数が素数かどうかを判定して、素数ならあとは全探索で$2^{N}$乗根を探す、といった感じでやりました。
さらに余談ですが、使うPは素数である必要はなく、Garnerで復元できる条件を満たすかつdigitNと互いに素であればなんでも良いです(多分)。

フェルマー小定理の参考
https://qiita.com/drken/items/6b4031ccbb2cab7436f3

3-3 負巡回畳み込みを使って0埋め領域を活用

ここまで触れて来ませんでしたが


    A=np.uint32([1,4,1,4,0,0,0,0])
    B=np.uint32([2,1,3,5,0,0,0,0])

の右半分の0は有効活用できないでしょうか。結論から言うとできます。ここを0以外にすると、結果の上位桁が一周回って下位桁に加算されます。
見てみたほうが早いので最初のfmt_smaple0.pyで以下の初期値にして試してみます。

    A=np.uint32([1,4,1,4,3,1,4,1])
    B=np.uint32([2,1,3,5,5,9,2,6])

結果

畳み込みの結果
C[ 0 ] = 97
C[ 1 ] = 78
C[ 2 ] = 73
C[ 3 ] = 87
C[ 4 ] = 61
C[ 5 ] = 77
C[ 6 ] = 87
C[ 7 ] = 67

これはこういうことになってます。
seijunkai.png

薄いピンクのところが1周回って畳み込まれています。これが正巡回畳み込みです。

次にピンクの上位桁の復元をみすえて負巡回畳み込みというのをします。同じA,Bを使ってFMTに使う$ω^{nk}$を$ω^{n(k+1/2)}$に変更します。すると

畳み込みの結果
C[ 0 ] = 164
C[ 1 ] = 197
C[ 2 ] = 202
C[ 3 ] = 222
C[ 4 ] = 15
C[ 5 ] = 25
C[ 6 ] = 75
C[ 7 ] = 67

hujunkai.png
こうなります。(画像とプログラムの出力は257剰余下でどちらも同じ値です。)

この正巡回と負巡回の結果を使えば
(正巡回-負巡回)/2 で上位桁
(正巡回+負巡回)/2 で下位桁
を復元できます。もちろん小数はありえない世界なので剰余下でやります。

少し戻って、さっきの$ω^{n(k+1/2)}$にするというところについてもう少し詳しく解説します。そもそも元のFMT関数はaを入力としてこのようなfを出力していました。

$f_k=a_0ω^{0}+a_1ω^{k}+a_2ω^{2k}....+a_{n-1}ω^{(n-1)k}$

負巡回ではこれをこうします。

$f_k=a_0ω^{0}+a_1ω^{k+(1/2)}+a_2ω^{2(k+(1/2))}....+a_{n-1}ω^{(n-1)(k+(1/2))}$

以下がこれをコードに書き起こしたものになります。
$ω=4$、法$P=257$として

#digitN=8
#MODP=257

#FMT関数と同じ。ナイーブO(N^2)実装
def naiveFMT(inputA):
    arrayA=np.zeros_like(inputA)
    for i in range(digitN):
        for j in range(digitN):
            arrayA[i]+=inputA[j]*ModExp(4,i*j)
            arrayA[i]%=MODP
    return arrayA[:]

#負巡回FMT
def naiveFMT_Neg(inputA):
    arrayA = np.zeros_like(inputA)
    for i in range(digitN):
        for j in range(digitN):
            arrayA[i]+=inputA[j]*ModExp(2,(i*2+1)*j)
            arrayA[i]%=MODP
    return arrayA[:]

こうなります。このやり方は$O(N^2)$の形ですが、実際の実装ではFMT関数部分を変えず、順変換の前にA[i],B[i]に$ω^{i/2}$をかけ、逆変換の後に$ω^{-i/2}$をかける処理を追加すれば良いです。この追加処理自体は線形時間なので負巡回の計算オーダーは正巡回と同じと言えます。
また今回は偶然$ω$の2乗根が存在したので大丈夫でしたが、$ω$の2乗根が存在しない場合は$ω^{1/2}$が計算できないので負巡回ができません。そのためにも原始$2^N$乗根の$2^N$がなるべく大きくなるものを探しておきたいものです。

負巡回畳み込みについて参考にさせて頂きました。
【FFTでの正巡回、負巡回畳み込みの話】
https://qiita.com/peria/items/cf4c4b72ebbeec7728af#fnref1
【FMT先駆者様のスライド】
http://ushiro.jp/research/fmt/sld021.htm

3-4 時間間引きFMTと周波数間引きFMTでビット逆順を相殺

時間間引きFMTとか周波数間引きFMTという言葉はありませんが、時間間引きFFTと周波数間引きFFTのFMT版ということです。FFTでもFMTでもバタフライ演算の前か後ろに配列ビット逆順の処理が必要です。

batah.png

batah2.png

稚拙なイメージ画像で申し訳ないですが、言いたいこととしては時間間引きFFTはバタフライ演算の前に配列ビット逆順が必要、周波数間引きFFTはバタフライ演算の後に配列ビット逆順が必要ということです。
上記ソースコードではFMT関数は時間間引きFFTを参考に書いてるので、バタフライ演算の前にビット逆順が必要になっています。逆変換も時間間引きです。ということは、順変換を周波数間引きにして、逆変換を時間間引きにすればビット逆順が相殺されます！途中にあるA[i]*B[i]のアダマール積は配列の順番関係ないので気にしなくていいです。

iranai.png

FMT：時間間引き順変換
iFMT：時間間引き逆変換
uFMT：周波数間引き順変換
iuFMT：周波数間引き逆変換
とし、負巡回も込みで実装すると

fmt_smaple4.py
import numpy as np

# aの逆元を返す mod m
def modinv(a,m):
    b = m
    u = 1
    v = 0
    while (b) :
        t = a // b
        a -= t * b
        tmp=a
        a=b
        b=tmp
        u -= t * v
        tmp = u
        u = v
        v = tmp
    u %= m
    if u < 0:
        u += m
    return u

#garnerのアルゴリズム、3変数ver
#64bit整数におさまらない可能性あり、pythonのint型を使用
def Garner(a_,b_,c_,ar_,br_,cr_):
    a = int(a_)
    b = int(b_)
    c = int(c_)
    ar = int(ar_)
    br = int(br_)
    cr = int(cr_)
    x = ar + (br - ar) * modinv(a,b) % b * a
    x = x + (cr - x % c) * modinv(a,c) % c * modinv(b,c) % c * a * b
    #x %= a*b*c #a<b<cのときこの処理は必要ない
    return x

class FMTClass():
    def __init__(self,MODP,MODP_WnSqrt):
        self.MODP=np.uint32(MODP)
        self.MODP_WnSqrt = np.uint32(MODP_WnSqrt)
        self.MODP_Wn = np.uint32(np.uint64(self.MODP_WnSqrt) * np.uint64(self.MODP_WnSqrt) % np.uint64(self.MODP))

    # MODP剰余下でaのb乗を返す
    # 追記1:np.uint32(pow(np.uint32(a),np.uint32(b),np.uint32(self.MODP)))を使ったほうが速い!
    def ModExp(self,a, b):
        ans = np.uint64(1)
        aa = np.uint64(a)
        while (b != 0):
            if (b % 2 == 1):
                ans = ans * aa % np.uint64(self.MODP)
            aa = aa * aa % np.uint64(self.MODP)
            b //= 2
        return np.uint32(ans)

    def FMT(self,inputA):
        arrayA=inputA.copy()
        for i in range(digit_level):
            loopCnt_Pow2 = 1 << i
            for idx in range(digitN // 2):
                t2 = idx % loopCnt_Pow2
                t0 = idx * 2 - t2
                t1 = t0 + loopCnt_Pow2
                arrayAt0 = arrayA[t0]
                arrayAt1 = arrayA[t1]
                w0 = self.ModExp(self.MODP_Wn, t2 * (1 << (digit_level - 1 - i)))
                w1 = np.uint32(np.uint64(arrayAt1) * np.uint64(w0) % np.uint64(self.MODP))
                arrayA[t1] = (arrayAt0 + self.MODP - w1) % self.MODP
                arrayA[t0] = (arrayAt0 + w1) % self.MODP
        return arrayA[:]

    def iFMT(self,inputA):
        arrayA = inputA.copy()
        for i in range(digit_level):
            loopCnt_Pow2 = 1 << i
            for idx in range(digitN // 2):
                t2 = idx % loopCnt_Pow2
                t0 = idx * 2 - t2
                t1 = t0 + loopCnt_Pow2
                arrayAt0 = arrayA[t0]
                arrayAt1 = arrayA[t1]
                w0 = self.ModExp(self.MODP_Wn, digitN - t2 * (1 << (digit_level - 1 - i)))
                w1 = np.uint32(np.uint64(arrayAt1) * np.uint64(w0) % np.uint64(self.MODP))
                arrayA[t1] = (arrayAt0 + self.MODP - w1) % self.MODP
                arrayA[t0] = (arrayAt0 + w1) % self.MODP
        return arrayA[:]

    def uFMT(self,inputA):
        arrayA = inputA.copy()
        for i in range(digit_level):
            loopCnt_Pow2 = 1 << (digit_level - i - 1)
            for idx in range(digitN // 2):
                t2 = idx % loopCnt_Pow2
                t0 = idx * 2 - t2
                t1 = t0 + loopCnt_Pow2
                arrayAt0 = arrayA[t0]
                arrayAt1 = arrayA[t1]
                w0 = self.ModExp(self.MODP_Wn, t2 * (1 << i))
                r0 = (arrayAt0 - arrayAt1 + self.MODP) % self.MODP
                r1 = (arrayAt0 + arrayAt1) % self.MODP
                w1 = np.uint32(np.uint64(r0) * np.uint64(w0) % np.uint64(self.MODP))
                arrayA[t1] = w1
                arrayA[t0] = r1
        return arrayA[:]

    def iuFMT(self,inputA):
        arrayA = inputA.copy()
        for i in range(digit_level):
            loopCnt_Pow2 = 1 << (digit_level - i - 1)
            for idx in range(digitN // 2):
                t2 = idx % loopCnt_Pow2
                t0 = idx * 2 - t2
                t1 = t0 + loopCnt_Pow2
                arrayAt0 = arrayA[t0]
                arrayAt1 = arrayA[t1]
                w0 = self.ModExp(self.MODP_Wn, digitN - t2 * (1 << i))
                r0 = (arrayAt0 - arrayAt1 + self.MODP) % self.MODP
                r1 = (arrayAt0 + arrayAt1) % self.MODP
                w1 = np.uint32(np.uint64(r0) * np.uint64(w0) % np.uint64(self.MODP))
                arrayA[t1] = w1
                arrayA[t0] = r1
        return arrayA[:]

    def DivN(self,inputA):
        digitN_r = modinv(digitN, int(self.MODP))  # digitNの逆元を計算
        ret = np.uint64(inputA[:]) * np.uint64(digitN_r) % np.uint64(self.MODP)
        return np.uint32(ret)

    def PreNegFMT(self,A):
        for i in range(digitN):
            A[i] = np.uint32(np.uint64(A[i]) * np.uint64(self.ModExp(self.MODP_WnSqrt, i)) % np.uint64(self.MODP))
        return A[:]

    def PostNegFMT(self,A):
        for i in range(digitN):
            A[i] = np.uint32(
                np.uint64(A[i]) * np.uint64(self.ModExp(self.MODP_WnSqrt, digitN * 2 - i)) % np.uint64(self.MODP))
        return A[:]

    #A*Bの畳み込み結果を返す
    #上位桁、下位桁を復元して2倍の要素数で返す
    def Convolution(self,A_,B_):
        #正巡回
        A = self.uFMT(A_%self.MODP)
        B = self.uFMT(B_%self.MODP)
        Cpos = np.uint64(A[:]) * np.uint64(B[:]) % np.uint64(self.MODP)  # 各要素どうし乗算
        Cpos = np.uint32(Cpos) #結果は32bit範囲内におさまってるので
        Cpos = self.iFMT(Cpos)  # 逆変換
        Cpos = self.DivN(Cpos)  # Nで割る

        # 負巡回
        A = self.PreNegFMT(A_ % self.MODP)#負巡回の前処理
        B = self.PreNegFMT(B_ % self.MODP)#負巡回の前処理
        A = self.uFMT(A)
        B = self.uFMT(B)
        Cneg = np.uint64(A[:]) * np.uint64(B[:]) % np.uint64(self.MODP)  # 各要素どうし乗算
        Cneg = np.uint32(Cneg)  # 結果は32bit範囲内におさまってるので
        Cneg = self.iFMT(Cneg)  # 逆変換
        Cneg = self.DivN(Cneg)  # Nで割る
        Cneg = self.PostNegFMT(Cneg)  # 負巡回の後処理

        #上位桁、下位桁 復元
        Elo = (Cpos + Cneg) % self.MODP
        Elo = (Elo + Elo % 2 * self.MODP) // 2

        Ehi = (Cpos - Cneg + self.MODP) % self.MODP
        Ehi = (Ehi + Ehi % 2 * self.MODP) // 2

        E = np.zeros(digitN*2).astype(np.uint32)
        E[:digitN] = Elo[:]
        E[digitN:] = Ehi[:]
        return E[:]

#検算用。A,Bの配列から一つの多倍長整数を作る。pythonのint型は乗算にKaratsuba法を使っているのでちょっと速い
def CalcAB(A,B):
    Alist = A.tolist()
    Blist = B.tolist()
    r=2**32
    for i in range(digit_level):
        for j in range(2**(digit_level-1-i)):
            Alist[j] = Alist[j * 2] + Alist[j * 2 + 1] * r
            Blist[j] = Blist[j * 2] + Blist[j * 2 + 1] * r
        r = r * r
    return Alist[0] * Blist[0]

#検算用。Eから一つの多倍長整数を作る
def CalcE(E):
    Elist = E.tolist()
    r = 2 ** 32
    for i in range(digit_level + 1):
        for j in range(2 ** (digit_level - i)):
            Elist[j] = Elist[j * 2] + Elist[j * 2 + 1] * r
        r = r * r
    return Elist[0]

def CreateWnSqrt(a,m,n_level):
    ans=a
    for i in range(n_level-digit_level-1):#負巡回で使うことを考え-1
        ans=ans*ans%m
    return ans

#E0,E1,E2を使ってgarnerでEを復元して繰り上がり処理
def Carrying(E0,E1,E2):
    E = np.zeros(digitN * 2).astype(np.uint32)
    upg = 0  # 繰り上がり分 python3のint型なのでオーバーフローの心配はない
    for i in range(digitN * 2):
        g = Garner(FMT_P0.MODP, FMT_P1.MODP, FMT_P2.MODP, E0[i], E1[i], E2[i])
        g += upg
        E[i] = np.uint32(g % 4294967296)
        upg = g//4294967296
    return E



if __name__=="__main__":
    digit_level = 16
    digitN = 1 << digit_level

    FMT_P0 = FMTClass(MODP=469762049 , MODP_WnSqrt=CreateWnSqrt(60733,469762049,26))
    FMT_P1 = FMTClass(MODP=1811939329, MODP_WnSqrt=CreateWnSqrt(59189,1811939329,26))
    FMT_P2 = FMTClass(MODP=2013265921, MODP_WnSqrt=CreateWnSqrt(52278,2013265921,27))

    # A,Bの初期化。前半後半ともにランダムな値をいれる。
    A_ = np.zeros(digitN).astype(np.uint32)
    A_[:] = np.uint32( np.random.randint(0, 2147483648, (digitN)) * 2
                                   + np.random.randint(0, 2, (digitN)))
    B_ = np.zeros(digitN).astype(np.uint32)
    B_[:] = np.uint32(np.random.randint(0, 2147483648, (digitN)) * 2
                                  + np.random.randint(0, 2, (digitN)))

    #p=469762049でやる
    E0 = FMT_P0.Convolution(A_,B_)

    #p=1811939329でやる
    E1 = FMT_P1.Convolution(A_, B_)

    #p=2013265921でやる
    E2 = FMT_P2.Convolution(A_, B_)

    #結果一部表示表示
    print("畳み込みの結果 mod","469762049",
          "\t\t畳み込みの結果 mod","1811939329",
          "\t\t畳み込みの結果 mod", "2013265921",
          "\t\tGarnerで復元")
    for i in range(10):
        print("E[",i,"] =",E0[i],
              "  \t\t\t\t", "E[", i, "] =", E1[i],
              " \t\t\t\t", "E[", i, "] =", E2[i],
              "\t\t\t\t",Garner(FMT_P0.MODP,FMT_P1.MODP,FMT_P2.MODP,E0[i],E1[i],E2[i])
              )
    #繰り上がり
    E=Carrying(E0,E1,E2)

    #検算
    print("検算開始")
    if CalcAB(A_, B_)==CalcE(E):
        print("正解！")
    else:
        print("ちがう！")

これでビット逆順の計算をなくし、A,Bの0埋め領域も消すことができました。また「FMT」関数と「iuFMT」関数は書いてはありますが実際には使われていません。

4 CUDAに移植

以上のコードをCUDAに移植します。私はPythonからCUDAをいじることが多いのでここではPyCUDAを使わせてください。

実行環境
GPU:RTX2080Ti
OS:Windows 10
Cuda compilation tools, release 10.1, V10.1.168
python 3.7.1 (Anaconda)
pyCUDA 2019.1 (CUDA 10.1に対応,pipコマンドでinstall)

(長いので折りたたんでいます。)

fmt_smaple5.py
fmt_smaple5.py
# FMT(高速剰余変換) PyCUDA版
# 3つの素数PのもとでFMTを行うサンプル
# garnerのアルゴリズムで最後に値を復元しているため、2の90.7乗まで正確に復元できる
# 検算関数もあり

# 追記1:これより1.2～1.8倍速いバージョンあり。(githubにある。剰余をモンゴメリ乗算に変えたやつ)
# 追記2:これより2～5倍速いバージョンあり。剰余をあらかじめ計算して参照するタイプ(sample6)
import pycuda.autoinit
import pycuda.driver as drv
from pycuda.compiler import SourceModule
import numpy as np
import time
import os

programid_garner=SourceModule("""
    #include "garner.cu"
""",include_dirs=[os.getcwd()])
kernel_garner=programid_garner.get_function("GarnerGPU")

class FMTClass():
    def __init__(self,MODP,MODP_WnSqrt):
        self.MODP=np.uint32(MODP)
        self.MODP_WnSqrt = np.uint32(MODP_WnSqrt)
        self.MODP_Wn = np.uint32(np.uint64(self.MODP_WnSqrt) * np.uint64(self.MODP_WnSqrt) % np.uint64(self.MODP))
        self.programid = SourceModule("""
            #define MODP (uint)(""" + str(MODP) + """)
                #include "fmt.cu"
            """, include_dirs=[os.getcwd()])
        self.kernel_FMT=self.programid.get_function("FMT")
        self.kernel_iFMT=self.programid.get_function("iFMT")
        self.kernel_uFMT=self.programid.get_function("uFMT")
        self.kernel_iuFMT=self.programid.get_function("iuFMT")
        self.kernel_Mul_i_i=self.programid.get_function("Mul_i_i")
        self.kernel_PostNegFMT=self.programid.get_function("PostNegFMT")
        self.kernel_PreNegFMT=self.programid.get_function("PreNegFMT")
        self.kernel_DivN=self.programid.get_function("DivN")
        self.kernel_PosNeg_To_HiLo=self.programid.get_function("PosNeg_To_HiLo")
        self.kernel_PostFMT_DivN_HiLo=self.programid.get_function("PostFMT_DivN_HiLo")


    def FMT(self,gpuMemA):
        for i in range(digit_level):
            self.kernel_FMT(gpuMemA, np.uint32(1 << i), self.MODP_Wn, np.uint32(1 << (digit_level - 1)),
                         grid=(gsz // lsz, 1, 1), block=(lsz, 1, 1))
        return

    def uFMT(self,gpuMemA):
        for i in range(digit_level):
            self.kernel_uFMT(gpuMemA, np.uint32((1<<(digit_level-1))>>i), self.MODP_Wn, np.uint32(1 << (digit_level - 1)),
                         grid=(gsz // lsz, 1, 1), block=(lsz, 1, 1))
        return

    def iFMT(self,gpuMemA):
        for i in range(digit_level):
            self.kernel_iFMT(gpuMemA, np.uint32(1 << i), self.MODP_Wn,
                         np.uint32(1 << (digit_level - 1)), grid=(gsz // lsz, 1, 1), block=(lsz, 1, 1))
        return

    def iuFMT(self,gpuMemA):
        for i in range(digit_level):
            self.kernel_iuFMT(gpuMemA, np.uint32(1<<(digit_level-1-i)), self.MODP_Wn,
                         np.uint32(1 << (digit_level - 1)), grid=(gsz // lsz, 1, 1), block=(lsz, 1, 1))
        return

    def Mul_i_i(self,gpuMemA,gpuMemB):
        self.kernel_Mul_i_i(gpuMemA, gpuMemB, grid=(gsz2 // lsz2, 1, 1), block=(lsz2, 1, 1))
        return

    def DivN(self,gpuMemA):
        self.kernel_DivN(gpuMemA, np.uint32(1<<digit_level), grid=(gsz2 // lsz2, 1, 1), block=(lsz2, 1, 1))
        return

    def PreNegFMT(self,gpuMemA,gpuMemB):
        self.kernel_PreNegFMT(gpuMemA,gpuMemB, self.MODP_WnSqrt,
                                    np.uint32(1<<digit_level), grid=(gsz2 // lsz2, 1, 1), block=(lsz2, 1, 1))
        return

    def PostNegFMT(self,gpuMemA):
        self.kernel_PostNegFMT(gpuMemA, self.MODP_WnSqrt,
                                     np.uint32(1<<digit_level), grid=(gsz2 // lsz2, 1, 1), block=(lsz2, 1, 1))
        return

    def PosNeg_To_HiLo(self,gpuMemE,gpuMemA,gpuMemB):
        self.kernel_PosNeg_To_HiLo(gpuMemE, gpuMemA,gpuMemB, np.uint32(1 << digit_level),
                                    grid=(gsz2 // lsz2, 1, 1), block=(lsz2, 1, 1))
        return

    def PostFMT_DivN_HiLo(self,gpuMemE,gpuMemA,gpuMemB):
        self.kernel_PostFMT_DivN_HiLo(gpuMemE, gpuMemA, gpuMemB, np.uint32(1 << digit_level),
                                 self.MODP_WnSqrt,grid=(gsz2 // lsz2, 1, 1), block=(lsz2, 1, 1))
        return

    # FMTで畳み込み乗算の結果を得る
    # A,Bは計算過程で内部が破壊されているので注意
    def Convolution(self,A,B):
        A_Neg = drv.mem_alloc(4 * digitN) #4=sizeof(uint)
        B_Neg = drv.mem_alloc(4 * digitN)  # 4=sizeof(uint)
        E = drv.mem_alloc(4 * digitN * 2)  # 4=sizeof(uint)

        self.PreNegFMT(A,A_Neg)#負順回用
        self.PreNegFMT(B,B_Neg)#負順回用
        self.uFMT(A)
        self.uFMT(B)
        self.uFMT(A_Neg)
        self.uFMT(B_Neg)
        self.Mul_i_i(A,B)#Bに結果が入る
        self.Mul_i_i(A_Neg,B_Neg)#Bに結果が入る
        self.iFMT(B)
        self.iFMT(B_Neg)
        self.PostNegFMT(B_Neg)
        self.DivN(B)#FFTみたくNで最後割らないといけない
        self.DivN(B_Neg)
        self.PosNeg_To_HiLo(E,B,B_Neg)#eに結果が入る
        #PostFMT_DivN_HiLo(E,B,B_Neg)   # 上4行のかわりにこれをつかってもよい

        drv.DeviceAllocation.free(A_Neg)
        drv.DeviceAllocation.free(B_Neg)
        return E


def Creategszlsz():
    gsz = 1 << (digit_level - 1)  # gpu global_work_size
    lsz = min(gsz, 256)  # gpu local_work_size
    gsz2 = digitN
    lsz2 = min(gsz2, 256)
    return gsz, lsz, gsz2, lsz2

#検算用。A,Bの配列から一つの多倍長整数を作る。pythonのint型は乗算にKaratsuba法を使っているのでちょっと速い
def CalcAB(A,B):
    Alist = A.tolist()
    Blist = B.tolist()
    r=2**32
    for i in range(digit_level):
        for j in range(2**(digit_level-1-i)):
            Alist[j] = Alist[j * 2] + Alist[j * 2 + 1] * r
            Blist[j] = Blist[j * 2] + Blist[j * 2 + 1] * r
        r = r * r
    return Alist[0] * Blist[0]

#検算用。Eから一つの多倍長整数を作る
def CalcE(E):
    Elist = E.tolist()
    r = 2 ** 32
    for i in range(digit_level + 1):
        for j in range(2 ** (digit_level - i)):
            Elist[j] = Elist[j * 2] + Elist[j * 2 + 1] * r
        r = r * r
    return Elist[0]

#検算
def AnswerCheck(host_A,host_B,host_E):
    print("検算開始")
    CPUres = CalcAB(host_A, host_B)
    GPUres = CalcE(host_E)
    if CPUres!=GPUres:
        print("not equal !!!!!!!")
        print(CPUres)
        print(GPUres)
    else:
        print("あってる")
    return

def InitializeAB():
    host_A = np.random.randint(2107483648, size=(digitN)).astype(np.uint32) * 2\
             + np.random.randint(2, size=(digitN)).astype(np.uint32)
    host_B = np.random.randint(2107483648, size=(digitN)).astype(np.uint32) * 2\
             + np.random.randint(2, size=(digitN)).astype(np.uint32)
    return host_A,host_B

def CreateWnSqrt(a,m,n_level):
    ans=a
    for i in range(n_level-digit_level-1):#負巡回で使うことを考え-1
        ans=ans*ans%m
    return ans

#E0,E1,E2を使ってgarnerでEを復元して繰り上がり処理。結果はGPU内
def Carrying(E0,E1,E2,E):
    gsz3 = 1 << (digit_level+1)
    lsz3 = min(gsz3, 256)
    kernel_garner(E0,E1,E2,E,np.uint32(digitN*2),
        grid=(gsz3 // lsz3, 1, 1),block=(lsz3, 1, 1))
    return





if __name__ == '__main__':
    digit_level = 25  #################################可変パラメーター############################# 最大25
    digitN = 1 << digit_level
    gsz, lsz, gsz2, lsz2 = Creategszlsz()#GPUカーネルのgridサイズblockサイズ計算

    print("A,Bの要素数=",digitN)
    FMT_P0 = FMTClass(MODP=469762049 , MODP_WnSqrt=CreateWnSqrt(60733,469762049,26))
    FMT_P1 = FMTClass(MODP=1811939329, MODP_WnSqrt=CreateWnSqrt(59189,1811939329,26))
    FMT_P2 = FMTClass(MODP=2013265921, MODP_WnSqrt=CreateWnSqrt(52278,2013265921,27))

    host_E = np.zeros(digitN * 2).astype(np.uint32)#結果格納用
    E = drv.to_device(host_E)  # gpuメモリ確保。かならず0に初期化しておかないといけない

    print("初期値生成")
    host_A, host_B = InitializeAB()
    A0 = drv.to_device(host_A) # gpuメモリ確保&転送
    B0 = drv.to_device(host_B) # gpuメモリ確保&転送
    A1 = drv.to_device(host_A)  # gpuメモリ確保&転送
    B1 = drv.to_device(host_B)  # gpuメモリ確保&転送
    A2 = drv.to_device(host_A)  # gpuメモリ確保&転送
    B2 = drv.to_device(host_B)  # gpuメモリ確保&転送

    print("GPU計算開始")
    stime = time.time()
    E0 = FMT_P0.Convolution(A0, B0)
    E1 = FMT_P1.Convolution(A1, B1)
    E2 = FMT_P2.Convolution(A2, B2)

    # ここからgarner&繰り上がり、E0,E1,E2から復元繰り上げしEに出力
    Carrying(E0, E1, E2 , E)
    drv.Context.synchronize() #タスク待ち命令
    print("GPU計算終了",time.time()-stime,"sec")

    # 結果をGPU→CPU
    drv.memcpy_dtoh(host_E, E)

    # 検算 ここは桁数が多いととても重い。digit_level = 20で約30秒
    #AnswerCheck(host_A, host_B, host_E)

fmt.cu
fmt.cu
#define ulong unsigned long long
#define uint unsigned int


//A*B%MODP
__device__ uint ABModC(uint a,uint b){
    ulong tmp=((ulong)(__umulhi(a,b)))*(1ULL<<32)+(ulong)(a*b);
    return (uint)(tmp%MODP);
}

//exp(a,b)%MODP
__device__ uint ModExp(uint a,uint b){
    ulong ans=1ULL;
    ulong aa=a;

    while(b!=0){
        if (b%2==1) ans=ans*aa%MODP;
        aa=aa*aa%MODP;
        b/=2;
    }
    return (uint)ans;
}


//逆変換後は、FFTでいうNで除算しないといけない。
// a/arrayLength mod P
__device__ uint DivN_f(uint a,uint arrayLength)
{
    uint as    =a/arrayLength;
    uint ar    =a-as*arrayLength;
    uint pn    =MODP/arrayLength;
    if (ar!=0){
        as+=(arrayLength-ar)*pn+1;
    }
    return as;
}




//arrayLength2 = arrayLength/2
__global__ void FMT(uint *arrayA,uint loopCnt_Pow2,uint omega,uint arrayLength2 ) {
    uint idx = threadIdx.x+blockIdx.x*256;
    uint t2 = idx%loopCnt_Pow2;
    uint t0 = idx*2-t2;
    uint t1 = t0+loopCnt_Pow2;
    uint w0;
    uint w1;
    uint arrayAt0=arrayA[t0];
    uint arrayAt1=arrayA[t1];
    uint r0;
    uint r1;
    w0=ModExp(omega,t2*(arrayLength2/loopCnt_Pow2));
    w1=ABModC(arrayAt1,w0);
    r0=arrayAt0-w1+MODP;
    r1=arrayAt0+w1;
    if (r0>=MODP){r0-=MODP;}
    if (r1>=MODP){r1-=MODP;}
    arrayA[t1]=r0;
    arrayA[t0]=r1;
}

__global__ void uFMT(uint *arrayA,uint loopCnt_Pow2,uint omega,uint arrayLength2 ) {
    uint idx = threadIdx.x+blockIdx.x*256;
    uint t2 = idx%loopCnt_Pow2;
    uint t0 = idx*2-t2;
    uint t1 = t0+loopCnt_Pow2;
    uint w0;
    uint w1;
    uint arrayAt0=arrayA[t0];
    uint arrayAt1=arrayA[t1];
    uint r0;
    uint r1;
    w0=ModExp(omega,t2*(arrayLength2/loopCnt_Pow2));
    r0=arrayAt0-arrayAt1+MODP;
    r1=arrayAt0+arrayAt1;
    if (r0>=MODP){r0-=MODP;}
    if (r1>=MODP){r1-=MODP;}
    w1=ABModC(r0,w0);
    arrayA[t1]=w1;
    arrayA[t0]=r1;  
}


__global__ void iFMT(uint *arrayA,uint loopCnt_Pow2,uint omega,uint arrayLength2 ) {
    uint idx = threadIdx.x+blockIdx.x*256;
    uint t2 = idx%loopCnt_Pow2;
    uint t0 = idx*2-t2;
    uint t1 = t0+loopCnt_Pow2;
    uint w0;
    uint w1;
    uint arrayAt0=arrayA[t0];
    uint arrayAt1=arrayA[t1];
    uint r0;
    uint r1;
    w0=ModExp(omega,arrayLength2*2-t2*(arrayLength2/loopCnt_Pow2));
    w1=ABModC(arrayAt1,w0);
    r0=arrayAt0-w1+MODP;
    r1=arrayAt0+w1;
    if (r0>=MODP){r0-=MODP;}
    if (r1>=MODP){r1-=MODP;}
    arrayA[t1]=r0;
    arrayA[t0]=r1;
}


__global__ void iuFMT(uint *arrayA,uint loopCnt_Pow2,uint omega,uint arrayLength2 ) {
    uint idx = threadIdx.x+blockIdx.x*256;
    uint t2 = idx%loopCnt_Pow2;
    uint t0 = idx*2-t2;
    uint t1 = t0+loopCnt_Pow2;
    uint w0;
    uint w1;
    uint arrayAt0=arrayA[t0];
    uint arrayAt1=arrayA[t1];
    uint r0;
    uint r1;
    w0=ModExp(omega,arrayLength2*2-t2*(arrayLength2/loopCnt_Pow2));
    r0=arrayAt0-arrayAt1+MODP;
    r1=arrayAt0+arrayAt1;
    if (r0>=MODP){r0-=MODP;}
    if (r1>=MODP){r1-=MODP;}
    w1=ABModC(r0,w0);
    arrayA[t1]=w1;
    arrayA[t0]=r1;
}


//同じ要素同士の掛け算
__global__ void Mul_i_i(uint *arrayA,uint *arrayB ) {
    uint idx = threadIdx.x+blockIdx.x*256;
    uint w0;
    w0=ABModC(arrayB[idx],arrayA[idx]);
    arrayB[idx]=w0;
}

//逆変換後のNで割るやつ。剰余下で割るには特殊処理が必要
__global__ void DivN(uint *arrayA,uint arrayLength ) {
    uint idx = threadIdx.x+blockIdx.x*256;
    arrayA[idx]=DivN_f(arrayA[idx],arrayLength);
}



//負巡回計算の前処理
//sqrt_omegaの2N乗が1 (mod P)
//a[0]*=ModExp(sqrt_omega,0)
//a[1]*=ModExp(sqrt_omega,1)
//a[2]*=ModExp(sqrt_omega,2)
//a[3]*=ModExp(sqrt_omega,3)
__global__ void PreNegFMT(uint *arrayA,uint *arrayB,uint sqrt_omega,uint arrayLength) {
    uint idx = threadIdx.x+blockIdx.x*256;
    arrayA[idx]%=MODP;
    uint w0=ModExp(sqrt_omega,idx);
    arrayB[idx]=ABModC(arrayA[idx],w0);
}

//負巡回計算の後処理
//sqrt_omegaの2N乗が1 (mod P)
//a[0]*=ModExp(sqrt_omega,-0)
//a[1]*=ModExp(sqrt_omega,-1)
//a[2]*=ModExp(sqrt_omega,-2)
//a[3]*=ModExp(sqrt_omega,-3)
__global__ void PostNegFMT(uint *arrayA,uint sqrt_omega,uint arrayLength) {
    uint idx = threadIdx.x+blockIdx.x*256;
    uint w0=ModExp(sqrt_omega,arrayLength*2-idx);
    arrayA[idx]=ABModC(arrayA[idx],w0);
}

//負巡回計算と正巡回計算結果から、上半分桁と下半分桁を求める
__global__ void PosNeg_To_HiLo(uint *arrayE,uint *arrayA,uint *arrayB,uint arrayLength) {
    uint idx = threadIdx.x+blockIdx.x*256;
    uint a=arrayA[idx];
    uint b=arrayB[idx];
    uint subab=(a-b+MODP);//まず(a-b)/2を求めたい
    uint flag=subab%2;
    subab-=MODP*((subab>=MODP)*2-1)*flag;//ここで絶対偶数になる
    subab/=2;//(a-b)/2 MOD Pを算出
    arrayE[idx+arrayLength]=subab;//上位桁は(a-b)/2 MOD P
    arrayE[idx]=a-subab+MODP*(a<subab);//a-((a-b)/2)=a/2+b/2 つまり(a+b)/2が下位桁
}


//vramへの書き込み回数を減らす目的に作った関数
//PostNegFMT関数とDivN関数とPosNeg_To_HiLo関数の統合版
__global__ void PostFMT_DivN_HiLo(uint *arrayE,uint *arrayA,uint *arrayB,uint arrayLength,uint sqrt_omega) {
    uint idx = threadIdx.x+blockIdx.x*256;
    uint a=arrayA[idx];
    uint b=arrayB[idx];

    //ここは負巡回の後処理計算部分
    uint w0=ModExp(sqrt_omega,idx+(idx%2)*arrayLength);
    b=ABModC(b,w0);

    //Nで除算する関数
    a=DivN_f(a,arrayLength);
    b=DivN_f(b,arrayLength);

    //あとは一緒
    uint subab=(a-b+MODP);//まず(a-b)/2を求めたい
    uint flag=subab%2;
    subab-=MODP*((subab>=MODP)*2-1)*flag;//ここで絶対偶数になる
    subab/=2;//(a-b)/2 MOD Pを算出
    arrayE[idx+arrayLength]=subab;//上位桁は(a-b)/2 MOD P
    arrayE[idx]=a-subab+MODP*(a<subab);//a-((a-b)/2)=a/2+b/2 つまり(a+b)/2が下位桁
}

garner.cu
garner.cu
#define ulong unsigned long long
#define uint unsigned int

#define MOD_P0 469762049LL
#define MOD_P1 1811939329LL
#define MOD_P2 2013265921LL


//３つの互いに素なPが与えられるので、それぞれの余りから元の値を復元したい
//このときPは全て固定なので剰余計算は全部決め打ちでいける

//E0～E2が入力、E3に出力
//繰り上がりも考慮
//arrayLength2=arrayE3の配列サイズ
__global__ void GarnerGPU(uint *arrayE0,uint *arrayE1,uint *arrayE2,uint *arrayE3,uint arrayLength2 ) {
    int idx = threadIdx.x+blockIdx.x*256;

    ulong ar=arrayE0[idx];
    ulong br=arrayE1[idx];
    ulong cr=arrayE2[idx];

    ulong x=ar;
    ulong brx=br-x+MOD_P1;
    if (brx>=MOD_P1)brx-=MOD_P1;
    x=x+(brx*1540148431)%MOD_P1*MOD_P0;
    //1540148431=modinv(MOD_P0,MOD_P1)
    //この時点でxはMOD_P1*MOD_P0以下であることが保証されている

    ulong crx=cr+MOD_P2-x%MOD_P2;
    if (crx>=MOD_P2)crx-=MOD_P2;
    ulong w1=(crx*1050399624)%MOD_P2;
    //1050399624=modinv(MOD_P0,MOD_P2) *modinv(MOD_P1,MOD_P2)%MOD_P2
    ulong w2=MOD_P0*MOD_P1;
    ulong out_lo=w1*w2;
    ulong out_hi=__umul64hi(w1,w2);

    if (out_lo>(out_lo+x)){
        out_hi++;
    }
    out_lo+=x;

    //ここから繰り上がり処理
    uint ui00_32=(uint)(out_lo%(1ULL<<32ULL));
    uint ui32_64=(uint)(out_lo/(1ULL<<32ULL));
    uint ui64_96=(uint)(out_hi%(1ULL<<32ULL));

    uint lastE3_0 = atomicAdd( &arrayE3[idx+0], ui00_32 );
    if ((lastE3_0+ui00_32)<lastE3_0){//繰り上がりを考慮
        ui32_64++;
        if (ui32_64==0)ui64_96++;
    }

    if (ui32_64!=0){
        uint lastE3_1 = atomicAdd( &arrayE3[idx+1], ui32_64 );
        if ((lastE3_1+ui32_64)<lastE3_1){//繰り上がりを考慮
            ui64_96++;//こいつがオーバーフローすることは絶対にない
        }
    }

    uint upflg=0;
    if (ui64_96!=0){
        uint lastE3_2 = atomicAdd( &arrayE3[idx+2], ui64_96 );
        if ((lastE3_2+ui64_96)<lastE3_2){//繰り上がりを考慮
            upflg++;
        }
    }

    uint lastE3_i;
    for(int i=idx+3;i<arrayLength2;i++){ //9999999+1みたいなとき用
        if (upflg==0)break;
        lastE3_i = atomicAdd( &arrayE3[i], upflg );
        if (lastE3_i==4294967295){
            upflg=1;
        }else{
            upflg=0;
        }
    }
}

3つのファイルは同じディレクトリに入れていれば動くはずです。
さてCUDAに移植する際にそのままは移植できないのでいくつか工夫する必要がありました。Garnerで復元するところはlong型を超える値が出てきてしまい(だいたい2進90桁)Python3 int型を使っていましたが、

    x = ar + (br - ar) * modinv(a,b) % b * a
    x = x + (cr - x % c) * modinv(a,c) % c * modinv(b,c) % c * a * b

このうちa,b,c,ar,br,crは31bitにおさまる数を使っているので、2行目の一番右のbを乗算するまでは確実に64bit整数の範囲内におさまります。ここでCUDAには__umul64hiという64bit整数×64bit整数の上位桁を計算する命令があるので、これを使うことで計算結果をulong型変数2つに余裕で格納することができます。

また、繰り上がり計算もGPUで並列計算しています。出力配列は1要素あたり32bitなので、Garnerで2進90桁を復元したら0-31桁、32-63桁、64-95桁と分割して各要素に加算します。

uint lastE3_1 = atomicAdd( &arrayE3[idx+1], ui32_64 );

ここではuint型のarrayE3[idx+1]にui32_64を加算したいという状況です。普通に加算処理を書くと、他のスレッドからもarrayE3[idx+1]にアクセスがあるのでデータ競合の問題が発生してしまいます。なのでここはAtomic演算を使って排他アクセスします。またatomicAddは加算前のarrayE3[idx+1]に入っている値を返してくれます。加算後の値というのは返って来ず、uint型でオーバーフローしたかもすぐには分かりません。ただ加算前の値+加算したい値はスレッド内でも計算できるので、arrayE3[idx+1]自体で繰り上がりが発生したかは判定できます(これはソースを見たほうが早いかもしれません)。そうしたらarrayE3[idx+2]に対してAtomic演算で加算処理しまたそこで繰り上がったらarrayE3[idx+3]に・・・といった処理をしていけば良いです。繰り上が続いていく確率は非常に低いので1スレッドあたり3,4回のAtomic演算で処理を終わらせることができます。Atomic演算自体はランダムアクセス程度の速度なので極端に遅くなる心配もないです。

また、CUDAでなくOpenCLでもmulhi命令やAtomic演算は扱えます。

さてそれでは速度比較をしたいと思います。

5 速度比較

5-1 Python vs CUDA

A,B=$2^N$桁(4294967296進数)
として
$C=A×B$
をCPU,GPUで計算して速度を測ります。
fmt_smaple4.py vs fmt_smaple5.pyをやっている感じです。
計測したところはA,Bのデータがデバイス(CPU,GPU)に確保された状態から、繰り上がり処理が完了したところまでです。

実行環境
CPU:core i9 9920X
RAM:DDR4 3600Mhz 16G*4 クアッドチャネル
GPU:RTX 2080Ti

処理時間(msec)

N Python PyCUDA
2 2.636432648 2.030849457
3 5.752801895 0.997066498
4 11.19756699 0.999450684
5 16.37482643 0.985145569
6 31.54706955 0.997781754
7 81.0508728 1.992464066
8 188.5592937 1.996278763
9 445.4622269 2.000808716
10 1047.091961 1.92117691
11 2431.799173 2.991437912
12 5611.189127 2.968549728
13 12762.6946 3.026247025
14 28899.42575 2.991199493
15 64907.24063 3.007650375
16 145894.0177 2.97832489
17 326845.052 4.023790359
18 732387.6169 8.923768997
19 1558503.907 16.954422
20 3391695.763 30.91740608
21 7350215.919 58.69293213
22 16059278.7 120.541811
23 34747664.01 222.6142883
24 75119683.56 473.2966423
25 161561500 940.9193993

speed_a.png
speed_b.png

な、なんということでしょう！最大で16万倍以上も高速化できました。$N=25$ではCUDAで0.94secに対しPythonで161000secかかっており歴然たる差です。GPGPU最高！！！
(定数倍高速化に目が行きがちですが、ちゃんとNが増えるにつき計算時間が2.3-4倍になっているところも注目です。)
正直なところここまで差がついたのは、私の書いたPythonコードが最適化されてないというのも大きいと思っています。無駄な型変換とか多い気がするし(追記1:コメント頂きpow関数使って2.6倍高速化できました)。もちろんCPU,GPUのハード性能の差も大きいはずです。インタプリタ言語特有の遅さみたいなのもありそうですし。なんだかいろんな要素が積み重なって16万倍の差につながっているようなので、まずはC++でも測ってみましょう。

5-2 Python vs C++ vs CUDA

C++ソースはこちら
https://github.com/toropippi/FMT_QiitaSample/tree/master/%E9%80%9F%E5%BA%A6%E6%AF%94%E8%BC%83%E7%94%A8C%2B%2B

処理時間(msec)

N Python C++ PyCUDA
2 2.6 20 2
3 5.8 8 1
4 11 7 1
5 16 24 1
6 32 9 1
7 81 9 2
8 189 9 2
9 445 21 2
10 1047 20 2
11 2432 33 3
12 5611 81 3
13 12763 144 3
14 28899 290 3
15 64907 656 3
16 145894 1456 3
17 326845 3255 4
18 732388 7183 8.9
19 1558504 15912 16.95
20 3391696 35121 30.92
21 7350216 77173 58.69
22 16059279 168893 120.5
23 34747664 367561 222.6
24 75119684 798197 473.3
25 161561500 1727122 940.9

speed_p_c_cu.png

Python vs C++で約90倍、C++ vs CUDAで約1800倍の差がつきました(C++はgcc -O3でコンパイル)。インタプリタ言語とコンパイラ言語の差もすごいです。こんなもんなんですかね？せいぜい10倍くらいの差かと思っていましたが。

ここで多倍長乗算の処理のうち何の処理が一番計算時間を占めているか気になったので計測してみました。C++での計算時間を参考にしています。

speed_circle.png

これを見ると計算時間のほとんどがFMT順変換と逆変換であることが分かります。いくら$O(NlogN)$のオーダーといえど二重ループになっているところであり、やはり重いですね。二重ループ内の処理はメモリアクセス(ランダムアクセス)と64bit整数の剰余計算がメインなのでそのどちらかが律速になっていると睨みました。そこでべき剰余だけを2倍量計算するコードでC++で計測してみたところ計算時間が約2倍になりました。つまりべき剰余の演算律速だったということです。

ボトルネックはおおむね特定できたので、先ほどのC++ vs CUDAの速度差についても考えてみたいと思います。現在1コア1スレッドSIMD化なしのC++プログラム vs CUDAという構図になっています。これではあまりにもCPUコード側が不利なので、CPU全コアフル稼働&SIMD化したらどこまで速くなるか気になります。
ところで剰余計算ってSIMD化でどこまで速くなるんですかね・・？除数が固定なら剰余が乗算に置き換えられてどうこうなると聞いたことがあるようなないような(追記1:剰余を乗算やビット演算に置き換えるのをモンゴメリ乗算というそうです。)

5-3 C++ vs OpenCL(on CPU)

CPUフル稼働できて私が使えるのはOpenCLしかないのでこれを使います。OpenCLはAMDのGPUでGPGPUやる時に登場しがちですが、CPU上でも動きます。私の経験ですがCPU上で動かしたときもかなり性能を引き出してくれる印象がありますので、性能比較で遊ぶにはもってこいです。
さっそく遊んでみましょう！

処理時間(msec)

N OpenCL(i9 9920X) C++
2 0 20
3 1 8
4 0 7
5 1 24
6 1 9
7 1 9
8 1 9
9 2 21
10 4 20
11 5 33
12 5 81
13 7 144
14 9 290
15 14 656
16 28 1456
17 55 3255
18 113 7183
19 240 15912
20 513 35121
21 1112 77173
22 2412 168893
23 5240 367561
24 11355 798197
25 24500 1727122

OpenCLはC++の70倍高速という結果です。i9 9920Xは12コアなので1コアあたりでも6-7倍高速化したことになります。これはSIMDの何がどう効いたのか良くわかりませんが・・・。(整数も早くなるんですかね)

5-4 OpenCL(on GPU) vs CUDA

せっかくなので同じ土俵でOpenCL vs CUDAの因縁の対決をやってみましょう。RTX2080Ti上でOpenCLとCUDAを動かしたときの速度はこうなりました。
処理時間(msec)

N OpenCL(RTX2080Ti) PyCUDA(RTX2080Ti)
2 1 2
3 0 1
4 1 1
5 1 1
6 1 1
7 1 2
8 1 2
9 1 2
10 1 2
11 1 3
12 1 3
13 2 3
14 1 3
15 2 3
16 3 3
17 4 4
18 6 9
19 13 17
20 26 31
21 54 59
22 113 121
23 242 223
24 435 473
25 930 941

ほんのわずかOpenCLのほうが速そうに見えますが。なんとなく引き分けでいいような気がします。
またOpenCLの結果からCPU(i9 9920X)とGPU(RTX2080Ti)の性能差が25倍くらいだということも分かりました。もちろんこの性能差というのは今回のプログラムに関してのみ言えることで、処理内容によっていくらでも変わりうるものです。

というわけでまとめると、16万倍の内訳は

Python vs C++ : コンパイラ言語の壁で90倍高速化
C++ vs OpenCL(on CPU) : マルチスレッド&SIMD化で70倍高速化
CPU vs GPU : ハードウェア性能差で25倍高速化

になります。

5-5 べき乗余のテーブル生成、参照で2～5倍高速化(2020/1追記分)

これまで二重ループ内のべき乗余$ω^{i}$を逐一計算してましたが、あらかじめ計算した配列を用意して参照するようにします。これでC++で5倍、PyCUDAで2倍速くなりました。
ソースはこちら
PyCUDA：https://github.com/toropippi/FMT_QiitaSample/blob/master/fmt_smaple6_qiita.py
OpenCL：https://github.com/toropippi/FMT_QiitaSample/tree/master/%E9%80%9F%E5%BA%A6%E6%AF%94%E8%BC%83%E7%94%A8OpenCL
C++：https://github.com/toropippi/FMT_QiitaSample/blob/master/%E9%80%9F%E5%BA%A6%E6%AF%94%E8%BC%83%E7%94%A8C%2B%2B/fmt_table.cpp

N C++テーブル参照なし C++テーブル参照あり
2 20 0
3 8 0
4 7 0
5 24 0
6 9 0
7 9 1
8 9 1
9 21 2
10 20 4
11 33 9
12 81 19
13 144 40
14 290 84
15 656 175
16 1456 369
17 3255 775
18 7183 1639
19 15912 3445
20 35121 7436
21 77173 15656
22 168893 32999
23 367561 72912
24 798197 174132
25 1727122 380685
N OpenCL(CPU)テーブル参照なし OpenCL(CPU)テーブル参照あり
2 0 0
3 1 0
4 0 1
5 1 1
6 1 1
7 1 1
8 1 1
9 2 1
10 4 2
11 5 3
12 5 3
13 7 3
14 9 4
15 14 5
16 28 7
17 55 12
18 113 19
19 240 35
20 513 74
21 1112 152
22 2412 323
23 5240 765
24 11355 2106
25 24500 4884
N PyCUDAテーブル参照なし PyCUDAテーブル参照あり
2 2 0
3 1 1
4 1 1
5 1 2
6 1 1
7 2 1
8 2 1
9 2 2
10 2 2
11 3 2
12 3 2
13 3 3
14 3 2
15 3 3
16 3 3
17 4 3
18 9 4
19 17 5
20 31 10
21 59 21
22 121 47
23 223 104
24 473 226
25 941 536

6 ガウス＝ルジャンドルのアルゴリズムで円周率を3億桁計算

いよいよ円周率を計算してみましょう。アルゴリズムはガウスルジャンドルwikiを使います。乗算以外に除算や平方根の計算がでてきますが、これはニュートン法を使えば良いです。
出力は16進数でtxtに出力とします。10進数は基数変換が大変そうなので次回以降実装とします。

ソースはこちら
https://github.com/toropippi/FMT_QiitaSample/blob/master/CalcPI_GPGPUfmt.py

speedamg.png

4294967296進数N桁 10進数N桁 計算時間(sec)
1024 9864.150898 1.105081319
2048 19728.3018 1.332427501
4096 39456.60359 1.591748238
8192 78913.20718 1.945914745
16384 157826.4144 2.207063675
32768 315652.8287 2.728562117
65536 631305.6575 3.273030281
131072 1262611.315 4.13736844
262144 2525222.63 5.635286331
524288 5050445.26 8.3940575
1048576 10100890.52 13.05427217
2097152 20201781.04 24.9399969
4194304 40403562.08 48.619494
8388608 80807124.16 98.77346897
16777216 161614248.3 190.9101679
33554432 323228496.6 423.4236162

だいたい円周率3億桁計算するのにGPUを使って420秒かかっています。y-cruncherだとCPUで10億桁200秒程度らしいのでまだまだですね・・。
アルゴリズム的にはy-cruncherが採用しているChudnovskyの公式が最速らしいのでガウスルジャンドルをそれに変えれば少しは速くなるでしょうか・・・？ただその場合binary splittingをどうGPUに最適化するかという問題になりそうです。というかなりました。それで挫折して今回ガウスルジャンドルを使った形です。
あとこのプログラムでは剰余計算が律速になっているようなのでそこをなんとか高速化したらもう少しy-cruncherに喰らいつけるのでは・・と思案してます。できるかもわかりませんが・・。
追記1:モンゴメリ乗算の愚直実装で1.2倍高速化できました。1.2倍・・・あまり速くならなかったけどこれ元々コンパイラが剰余をモンゴメリ乗算に展開してた可能性もある？
追記2:べき乗余テーブル参照にして2倍速くなりました(5-5に追記)

7 最後に

最後までお付き合い頂きありがとうございます。
実はFMTのGPU実装は6年前から挑戦していたのですが色んなポイントで挫折していました。6年も経てばQiita含め参考記事もいろいろ増えていて私自身、色んな面で助けられました。一方FMTで多倍長に挑戦した記事はあまりないように思い、同じことをやろうとしている人に少しでも役に立てばと思って記事にしてみました。単に多倍長の計算をするならGMPを使えばいいですが(それもそれで面白そうですが)、やはりGPUで計算するロマンと自分で実装する楽しさは何にも勝るため、何としてでも完成させてよかったです。参考にしてほしい反面、間違ったことを書いてるとそのまま広まってしまうのではという恐怖もあります。もし内容で気になったことがあればコメント頂ければ幸いです。
次回は未定です。

8 追加の参考文献とか

円周率jpのぺージ
http://xn--w6q13e505b.jp/

個人ブログ。C++実装コードあり
https://math314.hateblo.jp/entry/2015/05/07/014908

個人ブログ？C++実装コードあり
https://lumakernel.github.io/ecasdqina/math/FFT/NTT

個人ブログ。Python実装コードあり
https://smijake3.hatenablog.com/entry/2017/01/04/042136
https://tjkendev.github.io/procon-library/cpp/fft/fmt.html

ニュートン法で逆数と平方根を求める
http://www.finetune.co.jp/~lyuka/technote/fract/sqrt.html

モンゴメリ乗算
https://ikatakos.com/pot/programming_algorithm/number_theory/montgomery_multiplication

38
25
6

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
38
25