結果

問題 No.3020 ユークリッドの互除法・改
ユーザー iilj
提出日時 2025-02-14 21:41:35
言語 Python3
(3.13.1 + numpy 2.2.1 + scipy 1.14.1)
結果
AC  
実行時間 409 ms / 2,000 ms
コード長 5,302 bytes
コンパイル時間 155 ms
コンパイル使用メモリ 12,928 KB
実行使用メモリ 32,696 KB
最終ジャッジ日時 2025-02-14 21:41:48
合計ジャッジ時間 12,406 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge4
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 3
other AC * 21
権限があれば一括ダウンロードができます

ソースコード

diff # 

import numpy as np
from numpy.typing import NDArray


# https://qiita.com/yuji0001/items/64dc97cd4dcebf83d0a8


# 整数行列上の基本行列
def Eij(i: int, j: int, n: int) -> NDArray[np.int64]:
    E = np.eye(n, dtype=np.int64)
    E[i, i] = 0
    E[j, j] = 0
    E[i, j] = 1
    E[j, i] = 1
    return E


def Ei(i: int, n: int) -> NDArray[np.int64]:
    E = np.eye(n, dtype=np.int64)
    E[i, i] = -1
    return E


def Ec(i: int, j: int, c: int, n: int) -> NDArray[np.int64]:
    E = np.eye(n, dtype=np.int64)
    E[i, j] = c
    return E


# A[k:,k:]上の0以外の最小の絶対値をA[k,k]に移動する変換
def moveMN(
    A: NDArray[np.int64], k: int
) -> tuple[NDArray[np.int64], NDArray[np.int64], NDArray[np.int64]]:
    tmp_A = A[k:, k:]
    a = np.abs(tmp_A[tmp_A != 0]).min()
    i = np.where(np.abs(tmp_A) == a)[0][0] + k
    j = np.where(np.abs(tmp_A) == a)[1][0] + k
    P: NDArray[np.int64] = Eij(k, j, A.shape[1])
    invQ: NDArray[np.int64] = Eij(i, k, A.shape[0])
    B: NDArray[np.int64] = invQ.dot(A).dot(P)
    if B[k, k] < 0:
        Pi = Ei(k, A.shape[1])
        B = B.dot(Pi)
        P = P.dot(Pi)
    return invQ, B, P


# A[k,k]を使って，A[k+1:,k]を0に整える変換．(整いきれなかったら要素はA[k,k]より小さくなる)
def rowR(
    A: NDArray[np.int64], k: int
) -> tuple[NDArray[np.int64], NDArray[np.int64], NDArray[np.int64]]:
    B: NDArray[np.int64] = A.copy()
    invQ: NDArray[np.int64] = np.eye(A.shape[0], dtype=np.int64)
    P: NDArray[np.int64] = np.eye(A.shape[1], dtype=np.int64)
    for i in range(k + 1, A.shape[0]):
        q: int = A[i, k] // A[k, k]
        # 残渣
        # r: int = A[i, k] % A[k, k]
        invQi: NDArray[np.int64] = Ec(i, k, -q, A.shape[0])
        B = invQi.dot(B)
        invQ = invQi.dot(invQ)
    return invQ, B, P


# A[k,k]を使って，A[k,k+1]を0に整える変換．(整いきれなかったら要素はA[k,k]より小さくなる)
def colR(
    A: NDArray[np.int64], k: int
) -> tuple[NDArray[np.int64], NDArray[np.int64], NDArray[np.int64]]:
    B: NDArray[np.int64] = A.copy()
    invQ: NDArray[np.int64] = np.eye(A.shape[0], dtype=np.int64)
    P: NDArray[np.int64] = np.eye(A.shape[1], dtype=np.int64)
    for i in range(k + 1, A.shape[1]):
        q = A[k, i] // A[k, k]
        # 残渣
        # r = A[k, i] % A[k, k]
        Pi = Ec(k, i, -q, A.shape[1])
        B = B.dot(Pi)
        P = P.dot(Pi)
    return invQ, B, P


# A[k+1:,k+1:]においてA[k,k]で割り切れない要素A[i,j]をA[k,k]の残差に変換する変換
def remR(
    A: NDArray[np.int64], k: int
) -> tuple[NDArray[np.int64], NDArray[np.int64], NDArray[np.int64]]:
    invQ: NDArray[np.int64] = np.eye(A.shape[0], dtype=np.int64)
    P: NDArray[np.int64] = np.eye(A.shape[1], dtype=np.int64)
    #  Find i,j
    i = np.where(A[k + 1 :, k + 1 :] % A[k, k] != 0)[0][0] + k + 1
    j = np.where(A[k + 1 :, k + 1 :] % A[k, k] != 0)[1][0] + k + 1
    q = A[i, j] // A[k, k]
    # r = A[i, j] % A[k, k]
    invQi: NDArray[np.int64] = Ec(i, k, q, A.shape[0])
    Pi: NDArray[np.int64] = Ec(k, j, -1, A.shape[1])
    B: NDArray[np.int64] = invQi.dot(A).dot(Pi)
    P = P.dot(Pi)
    invQ = invQi.dot(invQ)
    return invQ, B, P


# Main Function
def Smith_Normalization(A: NDArray[np.int64]):
    invQ: NDArray[np.int64] = np.eye(A.shape[0], dtype=np.int64)
    P: NDArray[np.int64] = np.eye(A.shape[1], dtype=np.int64)
    A0 = A.copy()
    # limit of optimization
    N = 1000
    for k in range(min(A0.shape)):
        # If A0[k:,k:] is zero matrix, then stop calculation
        if np.sum(np.abs(A0[k:, k:])) == 0:
            break
        for i in range(N):
            if i == N - 1:
                print("Error: Time Out")
            # minimize A[k,k]
            invQi, A1, Pi = moveMN(A0, k)
            invQ = invQi.dot(invQ)
            P = P.dot(Pi)
            # make zero row A[k+1:,k]
            invQi, A2, Pi = rowR(A1, k)
            invQ = invQi.dot(invQ)
            P = P.dot(Pi)
            # if row A2[k+1:,k] is zero vector ?
            if np.abs(A2[k + 1 :, k]).sum() == 0:
                # make zero col A[k,k+1:]
                invQi, A3, Pi = colR(A2, k)
                invQ = invQi.dot(invQ)
                P = P.dot(Pi)
                # if col A3[k+1:,k] is zero vector ?
                if np.abs(A3[k, k + 1 :]).sum() == 0:
                    # A[k,k]|A[k+1:,k+1:]?
                    if np.sum(A3[k + 1 :, k + 1 :] % A3[k, k]) == 0:
                        A0 = A3.copy()
                        break
                    else:
                        # reduce A[k+1:,k+1:]
                        invQi, A0, Pi = remR(A3, k)
                        invQ = invQi.dot(invQ)
                        P = P.dot(Pi)
                else:
                    A0 = A3.copy()
            else:
                A0 = A2.copy()

    B = A0.copy().astype(int)
    P = P.astype(int)
    invQ = invQ.astype(int)
    return invQ, B, P


def main() -> None:
    A_raw = [list(map(int, input().split())) for _ in range(2)]
    # A_raw = [[100000000, 0], [0, 99999999]]
    # print(A_raw)
    # return

    # Check result
    A = np.array(A_raw, dtype=np.int64)
    _, B, _ = Smith_Normalization(A)
    # print(B)

    print(B[0, 0], B[1, 1])


if __name__ == "__main__":
    main()
0