AtCoder Beginner Contest 215 E

https://atcoder.jp/contests/abc215/tasks/abc215_e

ただのDPで、PyPy2だとあっさり通りました。
で、Python3で試したら、24ケース中13ケースでTLE。案の定Pythonでは遅いです。
手元でも試したら、実行時間は元の文字列の長さにだいたい比例するので、長さ少し変えて、

length Python2 PyPy2
1000 1.92s 0.25s
2000 3.83s 0.44s
4000 7.80s 0.84s

だいたいPyPy2のほうが9倍くらい速いですね。長いリストを使うとだいたいこんな感じです。
たいていの場合PyPyのほうが速いので、素直にPyPyを使いましょう。

AtCoder Regular Contest 118 D

https://atcoder.jp/contests/arc118/tasks/arc118_d

面白かったので解説してみます。

P = 13, a = 4, b = 5で考えてみます。
まず、グラフを描いてみましょう。

a倍していくと、
1 -> 4 -> 3 -> 12 -> 9 -> 10 -> 1
と周期6になります。
もう1つサークルがあります。
2 -> 8 -> 6 -> 11 -> 5 -> 7 -> 2

b倍していくと周期4が3つあります。
1 -> 5 -> 12 -> 8 -> 1
2 -> 10 -> 11 -> 3 -> 2
4 -> 7 -> 9 -> 6 -> 4

この2つの周期のサークルが絡み合っているので難しそうですよね。

さて、素数pに対して1 ~ p-1はaを適当に選ぶと周期p-1のサークルになります。p = 13に対してa = 6とすると、
1 -> 6 -> 10 -> 8 -> 9 -> 2 -> 12 -> 7 -> 3 -> 5 -> 4 -> 11 -> 1
周期12になります。
 g = 6と置くと、 4 \equiv g^{10}なので、a倍していくサークルは、
 1 \equiv g^0 \rightarrow g^{10} \rightarrow g^{20} \rightarrow g^{30} \rightarrow g^{40} \rightarrow g^{50} \rightarrow g^{60} \equiv 1
 2 \equiv g^5 \rightarrow g^{15} \rightarrow g^{25} \rightarrow g^{35} \rightarrow g^{45} \rightarrow g^{55} \rightarrow g^{65} \equiv 2

 5 \equiv g^9なので、b倍していくサークルは、
 1 \equiv g^0 \rightarrow g^9 \rightarrow g^{18} \rightarrow g^{27} \rightarrow g^{36} \equiv 1
 2 \equiv g^5 \rightarrow g^{14} \rightarrow g^{23} \rightarrow g^{32} \rightarrow g^{41} \equiv 2
 4 \equiv g^{10} \rightarrow g^{19} \rightarrow g^{28} \rightarrow g^{37} \rightarrow g^{46} \equiv 4

指数だけ見ましょう。そうすると、掛け算は足し算になって、12を法とする剰余を考えます。そして、指数を並べると、a倍を横、b倍を縦に並べて、

36 46 56 66 76 86 96
27 37 47 57 67 77 87
18 28 38 48 58 68 78
 9 19 29 39 49 59 69
 0 10 20 30 40 50 60

一番左と右、一番下と上は同じ数を表しています。つまり、上下左右に無限に伸びています。トーラスですね。
グラフの辺は、それぞれ数とその上下左右の数を結びます。
パスは、下の2行だけ考えて、

 9 → 19    29 → 39    49 → 59
 ↑    ↓    ↑     ↓    ↑     ↓
 0    10 → 20    30 → 40    50 → 60

とすればよいです。この方法は、横の周期が偶数だと成り立ちます。
ただし、b = 10だとパスはありません。なぜなら、 10 \equiv g^2なので、1にどうaとbを掛けても指数は偶数になってしまい、奇数は表せないからです。
 a \equiv g^{e_a}とすると、周期は p_a = (P-1)/\gcd{(P-1, e_a)}となります。逆に周期が分かっていたら、 \gcd{(P-1, e_a)} = (P-1)/p_aとなります。bの方の周期を p_bとすると、

 \gcd{(P-1, \gcd{(e_a, e_b)})} = \gcd{((P-1)/p_a, (P-1)/p_b)}

となります。これが1でないならパスは無しです。
周期はaを次々と掛けていけば得られると思いますが、周期はP-1の約数なので、約数の小さい順にべき乗の計算をすればよいでしょう。

最後に、横の周期を必ず偶数にできるかです。 p_a, p_bのどちらかが偶数になればいいので、どちらも奇数があり得るか考えます。
P-1はP=2以外では偶数です(P=2のときは、1 → 1というパスでOKです)。 p_a, p_bが共に奇数なら、 \gcd{((P-1)/p_a, (P-1)/p_b)}は偶数になります。

以下は、PyPy2のコードです。

# coding: utf-8
# Hamiltonian Cycle

from itertools import *
from fractions import gcd
import sys


#################### library ####################

def factorize_naive(n):
	def f(n, p0):
		if n <= 1:
			return ()
		
		for p in takewhile(lambda p: p * p <= n, count(p0)):
			if n % p == 0:
				return (p,) + f(n / p, p)
		else:
			return (n,)
	
	return list((k, len(list(v))) for k, v in groupby(f(n, 2)))

def divisors_f(fs):
	if len(fs) == 1:
		p, e = fs[0]
		for e1 in range(e+1):
			yield p**e1
	else:
		mid = len(fs) / 2
		ds1 = list(divisors_f(fs[:mid]))
		ds2 = list(divisors_f(fs[mid:]))
		for d1 in ds1:
			for d2 in ds2:
				yield d1 * d2

def period(e, p):
	fs = factorize_naive(p - 1)
	ds = sorted(divisors_f(fs))
	for d in ds:
		if pow(e, d, p) == 1:
			return d


#################### process ####################

def read_input():
	P, a, b = tuple(map(int, raw_input().split()))
	return (P, a, b)

def is_valid(path, P, a, b):
	if len(path) != P:
		return False
	
	s = set(path[:-1])
	if len(s) != P - 1:
		return False
	
	for v1, v2 in zip(path, path[1:]):
		if not (a * v1 % P == v2 or a * v2 % P == v1 or
				b * v1 % P == v2 or b * v2 % P == v1):
			return False
	else:
		return True

def Hamiltonian_Cycle(P, a, b):
	def generate_path(p1, e1, p2, e2):
		for i in range(p1):
			if i % 2 == 0:
				for j in range(p2):
					yield pow(e1, i, P) * pow(e2, j, P) % P
			else:
				for j in range(p2 - 1, -1, -1):
					yield pow(e1, i, P) * pow(e2, j, P) % P
		yield 1
	
	if P == 2:
		return [1, 1]
	
	pa = period(a, P)
	pb = period(b, P)
	if gcd((P - 1) / pa, (P - 1) / pb) != 1:
		return []
	
	if pa % 2 == 0:
		return list(generate_path(pa, a, (P - 1) / pa, b))
	else:
		return list(generate_path(pb, b, (P - 1) / pb, a))


#################### main ####################

P, a, b = read_input()
cycle = Hamiltonian_Cycle(P, a, b)
#print(is_valid(cycle, P, a, b))
if cycle:
	print('Yes')
	print(' '.join(map(str, cycle)))
else:
	print('No')

logの和の近似

\log{n!}の近似によくスターリングの公式が用いられます。
\log{n!} = n\log{n}-n
これで十分なことも多いですが、もう少しよい近似を考えてみましょう。

積分で近似して、
 \displaystyle S_1 = \int_2^{n+1}{\log{(x-\frac{1}{2}})dx}

この差を取ると、
 \displaystyle S - S_1 = \sum_{k=2}^n{(\log{k} - \int_k^{k+1}{\log{(x-\frac{1}{2}})dx})}

summationの中は、
 \log{k} - (k + \frac{1}{2})\log{(k + \frac{1}{2})} + (k - \frac{1}{2})\log{(k - \frac{1}{2})} + 1

 \log{(x+a)}テイラー展開すると、
 \log{(x+a)} = \log{x} + ax^{-1} - \frac{1}{2}a^2x^{-2} + \frac{1}{3}a^3x^{-3} + \cdots

 (x+a)\log{(x+a)}
 = (x+a)\log{x} + (x+a)(ax^{-1} - \frac{1}{2}a^2x^{-2} + \frac{1}{3}a^3x^{-3} + \cdots)
 = (x+a)\log{x} + a + a^2x^{-1} - \frac{1}{2}(a^2x^{-1}+a^3x^{-2}) + \frac{1}{3}(a^3x^{-2}+a^4x^{-3}) + \cdots
 = (x+a)\log{x} + a + \frac{1}{1\cdot 2}a^2x^{-1} - \frac{1}{2\cdot 3}a^3x^{-2} + \cdots

 (x+a)\log{(x+a)} - (x-a)\log{(x-a)}
 = 2a\log{x} + 2a - \frac{2}{2\cdot 3}a^3x^{-2} + \frac{2}{4\cdot 5}a^5x^{-4} + \cdots

だから、summationの中は、
 \frac{1}{2\cdot 3\cdot 2^2}k^{-2} + \frac{1}{4\cdot 5\cdot 2^4}k^{-4} + \cdots

 \displaystyle S - S_1 = \frac{1}{2\cdot 3\cdot 2^2}\sum_{k=2}^n{\frac{1}{k^2}} + \frac{1}{4\cdot 5\cdot 2^4}\sum_{k=2}^n{\frac{1}{k^4}} + \cdots

右辺第一項はリーマン・ゼータ関数

 \zeta(s) = \displaystyle \sum_{n=1}^{\infty}{\frac{1}{n^s}}

の特殊値

 \zeta(2) = \frac{\pi^2}{6}

が使えて、

 \displaystyle \sum_{k=2}^n{\frac{1}{k^2}} = \zeta(2) - 1 - \sum_{k=n+1}^{\infty}{\frac{1}{k^2}}

右辺第3項は上と同様な積分による近似が使えて、

 \displaystyle \sum_{k=n+1}^{\infty}{\frac{1}{k^2}} = \frac{1}{n+\frac{1}{2}} + O(n^{-3})

 S-S_1の右辺の第2項以降は、数値計算して、

from itertools import count

eps = 1e-22
v = []
for m in count(4, 2):
	a = 1.0 / (m*(m+1)*2**m)
	if a / 2**m < eps:
		break
	for k in count(2):
		b = a / k**m
		if b < eps:
			break
		v.append(b)

print(sum(sorted(v)))
0.000263942581577

 0.000263942581577 + O(n^{-3})

ということが分かりました。
以上をまとめると、

 S = (n+\frac{1}{2})\log{(n+\frac{1}{2})} - n + 1 - \frac{3}{2}\log{\frac{3}{2}} + \pi^2/144 + \frac{1}{24n+12} + 0.000263942581577 + O(n^{-3})

さらに定数項をまとめると、

 S = (n+\frac{1}{2})\log{(n+\frac{1}{2})} - n + \frac{1}{24n+12} + 0.46060519987133985 + O(n^{-3})

検証してみましょう。

from math import log

for n in (10, 20, 30):
	S = sum(log(m) for m in range(1, n+1))
	S2 = (n+0.5)*log(n+0.5)-n-1.0/(24*n+12)+0.41893853320467317
	print(n, S-S2)
>||

>||
10 2.0936225819667698e-06
20 2.8191458056880947e-07
30 8.56364152923561e-08

およそ-3乗で誤差が減っていっていることが分かると思います。

 n+\frac{1}{2}が気持ち悪かったら、そこも展開すると、

 S = n\log{n} - n + \frac{1}{2}\log{n} + \frac{1}{12n} + 0.91893853320467317 + O(n^{-3})

となりました。

新型コロナウイルスに関連した患者の死亡の統計

https://www.bousai.metro.tokyo.lg.jp/taisaku/saigai/1010035/1012790/index.html
ここから辿ると、死亡情報が例えば、

https://www.bousai.metro.tokyo.lg.jp/taisaku/saigai/1010035/1012790/1012818.html

のように得られます。年代・性別・居住地・診断日・死亡日が得られます。これに報告日をつけて、次のようなテーブルを得ました。

発表日 番号 年代 性別 居住地 診断日 死亡日
2021/1/14 1 90代 男性 都内 2020/12/9 2020/12/27
2021/1/14 2 70代 女性 都内 2020/12/27 2021/1/13
2021/1/14 3 80代 男性 都内 2020/12/31 2021/1/13

以下略

ただし、10/29以前はPDFでの情報提供のため、ここでは10/30以降の情報になっています。
これをRに読み込みます。

table <- read.csv2("death.csv", sep=",", colClasses=c("Date", NA, NA, NA, NA, "Date", "Date"))

発表日、診断日、死亡日のヒストグラムを描きます。

png("発表日.png")
hist(table[, 1], breaks="days", xlab="published day", freq=T, format="%m/%d")
dev.off()

f:id:inamori:20210116154127p:plain

png("診断日.png")
hist(table[, 6], breaks="days", xlab="diagnosis day", freq=T, format="%m/%d")
dev.off()

f:id:inamori:20210116154152p:plain

png("死亡日.png")
hist(table[, 7], breaks="days", xlab="death day", freq=T, format="%m/%d")
dev.off()

f:id:inamori:20210116154220p:plain

発表日より死亡日の方がバラツキが少ないことが分かります。
あと、死亡から発表のタイムラグと、診断から死亡までの日数が気になるので、これらもヒストグラムにします。

png("死亡から発表.png")
v <- as.integer(table[, 1] - table[, 7])
bins <- seq(min(v), max(v))
hist(v, breaks=bins, xlab="death to publish", freq=T)
dev.off()

f:id:inamori:20210116154658p:plain

png("診断から死亡.png")
v <- as.integer(table[, 7] - table[, 6])
bins <- seq(min(v), max(v))
hist(v, breaks=bins, xlab="diagnosis to death", freq=T)
dev.off()

f:id:inamori:20210116154717p:plain

診断から死亡はマイナスもけっこうあります。

PythonのライブラリをC++で作成する(5)

もう一つ、これはC/C++で行列の掛け算を書いたことがある人にはよく知られていると思いますが、乗数の方の行列を転置すると速くなります。
もう一度C++のコードを見てみましょう。

    Matrix  C(H, Vec(W));
    for(size_t i = 0U; i < H; ++i) {
        for(size_t j = 0U; j < W; ++j) {
            __int128    e = 0;
            for(size_t k = 0U; k < M; ++k)
                e += A[i][k] * B[k][j];
            C[i][j] = (long)(e % D);
        }
    }

一番内側のループでkを振っているので、行列Bを縦にスキャンしていくことになります。しかし、メモリは横に並んでいるので、飛び飛びに値を拾うことになり、キャッシュの中身が次々と入れ替わることになります。Bを転置すると、

    const Matrix    B_trans = transpose(B);
    Matrix  C(H, Vec(W));
    for(size_t i = 0U; i < H; ++i) {
        for(size_t j = 0U; j < W; ++j) {
            __int128    e = 0;
            for(size_t k = 0U; k < M; ++k)
                e += A[i][k] * B_trans[j][k];
            C[i][j] = (long)(e % D);
        }
    }

B_transは横にスキャンしています。これでキャッシュの入れ替わりが少なくなり、高速化されます。

実行してみましょう。

1134
256330968
12.208092451095581

これまでのバージョンの時間を並べてみましょう。

method time ratio
Python 1756.86s 225.8
PyPy 277.08s 22.7
C++ 206.85s 16.9
C++128 30.57s 2.5
C++transpose 12.21s 1

最後に、コードを表示します。

from distutils.core import setup, Extension

setup(name = 'mymatrix', version = '1.0.0', \
    ext_modules = [Extension('mymatrix', ['matrix.cpp'])])
// matrix.h
#include <Python.h>
#include <vector>

typedef std::vector<long>   Vec;
typedef std::vector<Vec>    Matrix;

static Vec make_vector(PyObject *obj_list);
static Matrix make_matrix(PyObject *args);
static PyObject *mul(PyObject *self, PyObject *args);
static Matrix mul_core(const Matrix& A, const Matrix& B, long D);
static Matrix mul_core128(const Matrix& A, const Matrix& B, long D);
static Matrix mul_core_trans(const Matrix& A, const Matrix& B, long D);
static Matrix transpose(const Matrix& A);
static PyObject *make_python_matrix(const Matrix& A);
static PyObject *make_python_vector(const Vec& v);
// matrix.cpp
#include <iostream>
#include "matrix.h"

using namespace std;

Vec make_vector(PyObject *obj_list) {
    if(!PyList_Check(obj_list))
        return Vec(0);
    
    Vec v;
    const int   M = PyList_Size(obj_list);
    for(int j = 0; j < M; ++j) {
        PyObject    *obj = PyList_GetItem(obj_list, (Py_ssize_t)j);
        if(!PyLong_Check(obj)) {
            printf("vector[%d] is not long.\n", j);
            return Vec(0);
        }
        
        const long  d = PyLong_AsLong(obj);
        v.push_back(d);
    }
    
    return v;
}

Matrix make_matrix(PyObject *obj_table) {
    if(!PyList_Check(obj_table))
        return Matrix(0);
    
    const int   L = PyList_Size(obj_table);
    Matrix  X(L);
    for(int i = 0; i < L; ++i) {
        PyObject    *obj = PyList_GetItem(obj_table, (Py_ssize_t)i);
        Vec v = make_vector(obj);
        if(v.empty())
            return Matrix(0);
        X[i] = v;
    }
    
    return X;
}

PyObject *mul(PyObject *self, PyObject *args) {
    PyObject    *obj1, *obj2;
    long    D;
    if(!PyArg_ParseTuple(args, "OOl", &obj1, &obj2, &D))
        return make_python_matrix(Matrix(0));
    
    const Matrix    E = make_matrix(obj2);
    const Matrix    A = make_matrix(obj1);
    const Matrix    B = make_matrix(obj2);
    Matrix  C = mul_core_trans(A, B, D);
    return make_python_matrix(C);
}

Matrix mul_core(const Matrix& A, const Matrix& B, long D) {
    const size_t    H = A.size();
    const size_t    M = B.size();
    const size_t    W = B.front().size();
    
    Matrix  C(H, Vec(W));
    for(size_t i = 0U; i < H; ++i) {
        for(size_t j = 0U; j < W; ++j) {
            for(size_t k = 0U; k < M; ++k)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % D;
        }
    }
    return C;
}

Matrix mul_core128(const Matrix& A, const Matrix& B, long D) {
    const size_t    H = A.size();
    const size_t    M = B.size();
    const size_t    W = B.front().size();
    
    Matrix  C(H, Vec(W));
    for(size_t i = 0U; i < H; ++i) {
        for(size_t j = 0U; j < W; ++j) {
            __int128    e = 0;
            for(size_t k = 0U; k < M; ++k)
                e += A[i][k] * B[k][j];
            C[i][j] = (long)(e % D);
        }
    }
    return C;
}

Matrix mul_core_trans(const Matrix& A, const Matrix& B, long D) {
    const size_t    H = A.size();
    const size_t    M = B.size();
    const size_t    W = B.front().size();
    
    const Matrix    B_trans = transpose(B);
    Matrix  C(H, Vec(W));
    for(size_t i = 0U; i < H; ++i) {
        for(size_t j = 0U; j < W; ++j) {
            __int128    e = 0;
            for(size_t k = 0U; k < M; ++k)
                e += A[i][k] * B_trans[j][k];
            C[i][j] = (long)(e % D);
        }
    }
    return C;
}

Matrix transpose(const Matrix& A) {
    const size_t    H = A.size();
    const size_t    W = A.front().size();
    Matrix  B(W, vector<long>(H));
    for(size_t i = 0U; i < W; ++i) {
        for(size_t j = 0U; j < H; ++j)
            B[i][j] = A[j][i];
    }
    return B;
}

PyObject *make_python_matrix(const Matrix& A) {
    const size_t    H = A.size();
    PyObject    *obj_mat = PyList_New(H);
    for(size_t i = 0U; i < H; ++i) {
        PyObject    *obj_list = make_python_vector(A[i]);
        PyList_SET_ITEM(obj_mat, i, obj_list);
    }
    return obj_mat;
}

PyObject *make_python_vector(const Vec& v) {
    const size_t    L = v.size();
    PyObject    *obj_list = PyList_New(L);
    for(size_t i = 0U; i < L; ++i)
        PyList_SET_ITEM(obj_list, i, PyLong_FromLong(v[i]));
    return obj_list;
}

static PyObject *print_matrix(PyObject *self, PyObject *args) {
    PyObject    *obj_table;
    if(!PyArg_ParseTuple(args, "O", &obj_table))
        return make_python_matrix(Matrix(0));
    
    Matrix  X = make_matrix(obj_table);
    for(auto p = X.begin(); p != X.end(); ++p) {
        for(auto q = p->begin(); q != p->end(); ++q)
            cout << *q << ' ';
        cout << endl;
    }
    
    return PyLong_FromLong(0);
}

static PyMethodDef mymatrixMethods[] = {
    { "print_matrix", print_matrix, METH_VARARGS },
    { "mul", mul, METH_VARARGS },
    { NULL }
};

static struct PyModuleDef mymatrix = {
    PyModuleDef_HEAD_INIT, "mymatrix", "Python3 C API Module",
    -1, mymatrixMethods
};

PyMODINIT_FUNC PyInit_mymatrix() {
    return PyModule_Create(&mymatrix);
}
# panel.py
# coding: utf-8
from itertools import *
from collections import Counter
import sys
import time
import mymatrix


#################### Matrix ####################

def mat_pow(M, e, D):
    if e == 1:
        return M
    elif e % 2 == 1:
        return mymatrix.mul(M, mat_pow(M, e-1, D), D)
    else:
        A = mat_pow(M, e//2, D)
        return mymatrix.mul(A, A, D)


#################### process ####################

def f(H, W, D):
    def normalize(diagram):
        min_l = min(diagram)
        return tuple(row - min_l for row in diagram)
    
    def nexts(diagram):
        y = next(y for y, row in enumerate(diagram) if row == 0)
        yield diagram[:y] + (2,) + diagram[y+1:]
        if y+1 < H and diagram[y+1] == 0:
            yield diagram[:y] + (1, 1) + diagram[y+2:]
    
    def walk(diagram):
        for diagram1 in nexts(diagram):
            diagram2 = normalize(diagram1)
            if diagram2 not in set_diagram:
                set_diagram.add(diagram2)
                walk(diagram2)
    
    def make_graph(H):
        graph = Counter()       # { (orig, dest): counter }
        for diagram in set_diagram:
            for diagram1 in nexts(diagram):
                graph[(diagram, normalize(diagram1))] += 1
        return graph
    
    def make_matrix(H):
        graph = make_graph(H)
        M = [ [ graph[(orig, dest)] for orig in set_diagram ]
                                    for dest in set_diagram ]
        row_to_index = dict((diagram, i) for i, diagram in enumerate(set_diagram))
        return (M, row_to_index)
    
    init_diagram = (0,)*H
    set_diagram = set([init_diagram])
    walk(init_diagram)
    M, row_to_index = make_matrix(H)
    print(len(M))
    M_pow = mat_pow(M, W*H//2, D)
    i = row_to_index[init_diagram]
    return M_pow[i][i]


#################### main ####################

t0 = time.time()
D = 10**9+7
H = int(sys.argv[1])
W = int(sys.argv[2])
print(f(H, W, D))
print(time.time() - t0)

PythonのライブラリをC++で作成する(4)

なぜC++で書いたライブラリは遅いのでしょう。
Pythonのコードは、

    for i in range(H):
        for j in range(W):
            C[i][j] = int(sum(A[i][k] * B[k][j] for k in range(M)) % D)
    return C

C++のコードは、

    for(size_t i = 0U; i < H; ++i) {
        for(size_t j = 0U; j < W; ++j) {
            for(size_t k = 0U; k < M; ++k)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % D;
        }
    }

Pythonのコードは総和を取ったあとにDの剰余を取っていますが、C++のコードは掛け算して足し算するたびにDの剰余を取っています。
なぜこうなっているかというと、競技プログラミングではDが32ビットに近い数が選ばれるため、すぐに64ビットではオーバーフローしてしまうからです。一方、Pythonでは多倍長整数に移行するため、オーバーフローを考慮しなくてよく、また多倍長整数の加算は64ビット整数の除算よりずっと速いです。
C++では64ビットより長い整数を使えばよいですが、幸いにしてGCCには__int128という128ビット整数があり、これを利用すれば速くなります。

PyObject *mul(PyObject *self, PyObject *args) {
    ...
    Matrix  C = mul_core128(A, B, D);
    return make_python_matrix(C);
}

Matrix mul_core128(const Matrix& A, const Matrix& B, long D) {
    const size_t    H = A.size();
    const size_t    M = B.size();
    const size_t    W = B.front().size();
    
    Matrix  C(H, Vec(W));
    for(size_t i = 0U; i < H; ++i) {
        for(size_t j = 0U; j < W; ++j) {
            __int128    e = 0;
            for(size_t k = 0U; k < M; ++k)
                e += A[i][k] * B[k][j];
            C[i][j] = (long)(e % D);
        }
    }
    return C;
}

実行してみましょう。

1134
256330968
30.565128803253174

変更前に対して6.8倍速くなりました。
速度アップの方法がもう一つあります。