C++で素数Pの倍数となる等差数列の最初の項を求める方法

概要

初項 A と公差 d からなる等差数列(AP)と、ある素数 P が与えられたとき、この等差数列の中で素数 P の倍数となる最初の項の位置(インデックス)を求めるのが本記事の目的です。

入力

A = 3, d = 4, P = 5

出力

3

説明

この等差数列の第4項が、素数 5 の倍数になっています。

• 第1項 = 3
• 第2項 = 3 + 4 = 7
• 第3項 = 3 + 2×4 = 11
• 第4項 = 3 + 3×4 = 15

15 = 5 × 3 であるため、第4項(インデックス 3)が最初の P の倍数となります。

解法の考え方

N番目の項を A_N とすると、次のように表せます。

A_N = A + (N−1)×d

A_N が P の倍数であるという条件は、次の式で書けます。

A + (N−1)×d = l×P(ここで l は定数)

ここで、A を A % P、d を d % P に置き換えると、(N−1)×d = l×P − A となります。

右辺に P を加算・減算すると、

(N−1)×d = P(l−1) + (P − A)

このとき、A は P より小さい値(A % P)に置き換えられているため、P − A は非負の数になります。最後に両辺の mod を取ると、

((N−1)×d) % P = (P − A) % P、すなわち ((N−1)×d) % P = P − A

ここで、(d×Y) % P = 1 を満たす P 未満の Y を考えます。この Y は、P を法とする d のモジュラ逆数(逆元)と呼ばれます。モジュラ逆数はフェルマーの小定理より、d^P−2 % P として高速に計算できます。

したがって、求める答え N は次の式で表されます。

N = ((Y × (P − A)) % P) + 1

実装例(C++)

#include <bits/stdc++.h>
using namespace std;
// (x1^y1)%p1 を O(log y1) で計算する反復関数
int power(int x1, int y1, int p1){
    // 結果の初期化
    int res1 = 1;
    // x1 が p 以上の場合は更新する
    x1 = x1 % p1;
    while (y1 > 0) {
        // y1 が奇数の場合、x1 を結果に乗算する
        if (y1 & 1)
            res1 = (res1 * x1) % p1;
        // この時点で y1 は必ず偶数
        y1 = y1 >> 1; // y1 = y1/2
        x1 = (x1 * x1) % p1;
    }
    return res1;
}
// P の倍数となる最初の項の位置を求める関数
int NearestElement1(int A, int d, int P){
    // 基本条件
    if (A == 0)
        return 0;
    else if (d == 0)
        return -1;
    else {
        int Y = power(d, P - 2, P);
        return (Y * (P - A)) % P;
    }
}
// ドライバーコード
int main(){
    int A = 3, d = 4, P = 5;
    // A と d の両方を mod する
    A %= P;
    d %= P;
    // 関数呼び出し
    cout << NearestElement1(A, d, P);
    return 0;
}

出力

3

計算量

モジュラ冪乗の計算に O(log P) かかるため、全体の時間計算量は O(log P) となります。これは項を順番に1つずつ確認する O(P) の線形探索よりもはるかに効率的です。

まとめ

等差数列と素数が与えられた場合、モジュラ逆数を利用することで、P の倍数となる最初の項の位置を対数時間で求められます。初項がすでに 0 の場合は 0 を、公差が 0 の場合は解が存在しないため -1 を返すという境界条件の処理も実装のポイントです。