前言：刷题时意外碰到这题，初见时以为是简单的数论题，求出通项公式再分析即可，但实际又由于出现根号问题无法直接分析。求助好友后才发现是图论题，觉得这种建图方式很有趣，故记录下来，并延申。

题意

ABC446E

题意就是对于一个序列有：

求对于集合 $0\le x \ , \ y \le M-1$ ，x和y不同取值下，无限长的序列中不存在 $M \mid s_i$ （即 $s_i$ 不被 M 整除）的(x,y)取值个数。

思路

首先这个通项公式很好写出，就是：

其中 $\lambda_1 \ , \ \lambda_2$ 是方程 $x^2-Ax-B=0$ 的根， $c_1 \ , \ c_2$ 是由 $s_1=x$ 和 $s_2=y$ 决定的常数。

但是很显然是这个底数 $\lambda$ 不一定是整数，比如著名的斐波那契数列，其底数就是带有根号的黄金分割数，直接分析还需要考虑用二项定理展开，然后根据具体的 $c_1$ 和 $c_2$ 的值，消去根号才能分析，实属复杂。

注意到这里的模数M只有1000，尝试暴力会发现，我们需要分别枚举 $x$ 和 $y$ ，然后对于每个独立的(x,y)，计算在前M*M项数中是否存在%M==0

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
for(int i=0;i<M;i++)for(int j=0;j<M;j++)
{
    bool ok=true;
    vector<int> s(M*M+1);
    s[1]=i, s[2]=j;
    for(int k=3; k<=M*M; k++)
    {
        s[k]=((ll)A*s[k-1]+(ll)B*s[k-2])%M;
        if(s[k]==0)
        {  
            ok=false;
            break;
        }
    }
    if(ok) ans++;
}

但是这样的复杂度是 $O(M^4)$ ，显然TLE。

图

这里考虑所有的x和y的取值后，序列 ${s}$ 的值必然在模数 $M$ 内
我们作状态 $(u,v):=(s[i-1],s[i])$ ，那么有转移 $(s{i-2},s{i-1}) \to (s_{i-1},s_i)$ ，也就是每个节点有且仅有一条出边（不一定只有一条入边），即外向基环树

那么不同对的(x,y)又应该如何考虑呢？

我们注意到对于 A=1, B=2;
我们取 $(x=0,y=1)$ 有

注意这个节点 $(1,1)$ ，如果我们取 $(x=1,y=1)$ 在数值上也可以得到这个结果 $(s[1]=1,[2]=1)$
也就是说，对于该空间 $M^2$ 内，节点的数值既可以表示 $(x,y)$ 对，同时也是 $(s[i-1],s[i])$

对于该空间 $M^2$ 具体的节点是：

当然其具体的连接关系为:

题目求有多少个(x,y)产生的序列不含M的倍数，在这个状态空间中也就是确定一个起点(x,y)沿着其出边一直走， 并使其不能碰到x=0或者y=0，统计这样的路径。但是这样直接遍历每个节点（上面的方法），复杂度还是 $O(M^4)$ ，显然会超。

这里说明下一些应该知道的关系

  • $M^2=$ 经过含0节点的路径数 $+$ 没经过含0节点的路径数
  • 对于所有的 $(x,y)$ 或 $(x,0)$ 或 $(0,0)$ 作为起点的序列必然会经过含0节点
  • 如果一个节点的后继是含0节点 那么其前驱的后继的后继也是会经过含0节点

那么必然地，

  • 对于从0节点出发，反着走，经过的每个节点所对应的序列，都会经过0节点
    • 也就是从所有的0节点反着出发，经过的每个节点都视为必然有一条路径经过
      • 即经过每个节点计数+1

我们只要 $M^2$ 地遍历节点构造反图，然后从含0节点bfs搜索计数标记节点，然后 $M^2-\text{cnt}$ 即可

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
inline int encode(int x, int y, int m)
{
	return x*m+y;
}

int main()
{
	ios::sync_with_stdio(false);
	cin.tie(0);
	cout.tie(0);
	int m,a,b; cin >> m >> a >> b;
	// 这里使用编码存(x,y) x*m+y
	vector<vector<int> > rg(m*m);
	for(int i=0; i<m; i++)
	{
		for(int j=0; j<m; j++)
		{
			int x1=i, y1=j;
			int x2=j, y2=((ll)a*j+(ll)b*i)%m;
			int u=encode(x1,y1,m);
			int v=encode(x2,y2,m);
			rg[v].push_back(u);
		}
	}
	vector<bool> vis(m*m,false);
	queue<int> q;
	
	for(int k=0; k<m; k++)
	{
		int u1=encode(0,k,m);
		if(!vis[u1])
		{
			vis[u1]=true;
			q.push(u1);
		}
		int u2=encode(k,0,m);
		if(!vis[u2])
		{
			vis[u2]=true;
			q.push(u2);
		}
	}
	
	int cnt=0;
	while(!q.empty())
	{
		int u=q.front(); q.pop();
		cnt++;
		for(int v:rg[u])if(!vis[v])
		{
			vis[v]=true;
			q.push(v);
		}
	}
	
	cout << m*m-cnt << endl;
	
	return 0;
}

总结

本题很巧妙地使用了状态空间将序列的集合和(x,y)初始选点的集合合并起来，使得在搜索图（序列结果）中恰好也统计了选点的个数，很妙！

而且注意到每个状态到下一个状态是唯一的，即每个结点只有一条出边，整个图由若干个 基环树 Functional Components 组成，这意味着从任何点出发，最终都会进入一个环。

扩展

下面引出一些数论问题和图论问题的重合部分：

序列的周期与相遇问题

给定 $x_{n+1} := (Ax_n + C) \pmod M$，求序列从何时开始进入循环，循环长度是多少？

数论解法

该问题如果使用数论直接解的话十分复杂，这里只给出大概解法：

1) $C=0$ 时

如果 $C=0$ 即讨论 $x_n = A^n x_0 \pmod M$ 的周期和进入周期情况。

我们先把 $x_0$ 消掉，考虑同余式 $A^nx_0 \equiv A^{n+k}x_0 \pmod M$ ，那么同时除以其最小公倍数 $g=\gcd(x_0,M)$ 得到 ：

注意到 $\gcd(\frac{x_0}{g} ,\frac{M}{g})=1$ ，所以 $\frac{x_0}{g}$ 在模 $\frac{M}{g}$ 下有逆元，上式同时乘以 $\frac{x_0}{g}$ 的逆元，得到：

然后考虑模m意义下对A幂的影响，提取出关键因子，考虑 $m$ 的质因数分解，设 $m=p_1^{q_1}p_2^{q_2}\cdots p_n^{q_n}$ ，并定：

也就是将 $M$ 分解为含整除 $A$ 的质因数和不含整除 $A$ 的质因数，那么 $(*)$ 式成立当且仅当（CRT）：

  • 对于模 $m_2$ 的项， $\gcd(A,m_2)=1$ ，由欧拉定理有：

    即在对于这个模数来说，序列 $A \pmod{m_2}$ 是循环的，且最小正周期 $t$ 为 $A \pmod{m_2}$ 的阶，记作：

    这里阶满足 $\delta{m_2}(A) | \varphi(m_2)$ 即阶为 $m_2$ 的欧拉函数的约数，这样可以用程序快速求出。
    例如先初始化阶等于其欧拉函数，尝试是否满足 $A^{\delta    {m_2}(A)} \equiv 1 \pmod {m_2}$ ，然后不断除以 $\varphi(m_2)$ 的每个质因子，直到不能满足余数为1，此时上个阶即为所求。

  • 对于模 $m_1$ 的项，由于其所有素因子均能整除 $A$，设存在足够大的 $n$ 使得 $A^n$ 包含 $m_1$ 所需的所有素因子幂次。
    即此时 $n$ 满足：

    那么对于之后的 $n+k \ , \ k>1$ ， $A^{n+k} \pmod{m_1}$ 仍然满足上述条件。
    我们找到其预周期也就是找到最小的 $t$，使得对于所有 $n \ge t$，都有：

    我们记 $v_p(x)$ 为 $x$ 的 $p$ 进赋值（即x中质因子p的最高幂次），则对于 $A^n \equiv 0 \pmod{p^{v_p(m_1)}}$ 满足：

    因此，进入周期 $T$ 的最小预周期 $t$ 为：

    若 $m_1=1$ （即 $\gcd(A,m)=1$ ），则 $t=0$ ，这也印证了当A和m互质时，直接进入大循环。

2) $C \ne 0$ 时

引理：线性同余方程的齐次化 Homogenization Lemma

给定线性同余递归序列 $x_{n+1} \equiv Ax_n + C \pmod M$ ，若满足 $\gcd(A-1, M) = 1$ ，则该序列可以通过坐标平移 $y_n = x_n - x^*$ 完美等价于齐次序列（即 $C=0$ 的形式）：

其中 $x^*$ 为递归映射的不动点。

证明:

我们需要找到一个常数 $x^$ ，使得当 $x_n = x^$ 时， $x_{n+1}$ 也等于 $x^*$ 。根据递推式：

整理得：

又 $\gcd(1-A, M) =\gcd(A-1, M) = 1$ ，那么 $(1-A)$ 在模M下有逆元，即：

那么将坐标平移 $y_n = x_n - x^*$ ，并带入原式，得：

也就是只要找到 $x^*$ ，就可以将原递推式等价于齐次递推式，从而利用上述方法求解。

由引理，可得如果 $\gcd(A-1, M) = 1$ 那么我们直接将递推式齐次化（也就是使得C=0），然后化归为前面的问题，求解即可。

但是如果 $\gcd(A-1, M) \ne 1$ ，即 $(A-1)$ 没有逆元，我们可以尝试扩大模数，令 $zn:=(A-1)x_n+C$ ，那么 $z{n+1}$ 可以表示为：

神奇地，$z_{n}$ 的递推关系竟然是齐次的，由于z相对于x被扩大了(A-1)，那么其模数也应该扩大到 $M(A-1)$ :

这与前面齐次的解法一样。

图论解法

这个问题就很显然了，我们直接作状态空间 $xn$ 由于是模M的，所以 $x_n={ 0,1,2\cdots M-1 }$ ，其连接关系就是 $x{n+1} = (Ax_n + C) \pmod M$ 。
我们要求的是进入周期的步长 $t$ ，即从 $x_0=C \pmod M$ 出发，第一次到达环的步长。然后环的长度就是周期 $T$ 。

程序
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
using ll = long long;

ll A,C,M;

inline ll f(ll x)
{
	return (A*x+C)%M;
}

signed main()
{
	cin >> A >> C >> M;
	// floyd
	ll sl=f(0);
	ll fa=f(sx); // 快慢指针
	while(sl!=fa)
	{
		sl=f(sl);
		fa=f(f(fa));
	}
	// sx,nx都在环上相遇
	// 且i,j到环的入口距离相同
	ll t=0; // 预周期
	ll i=f(0), j=sl;
	// 让i处于x0的位置，j在环入口，不断移动直到相遇
	while(i!=j)
	{
		i=f(i);
		j=f(j);
		t++;
	}
	ll T=1; // 周期
	j=f(i);  // j从环入口下一个结点开始
	while(i!=j) j=f(j), T++;
	
	cout << "pre period:" << t << endl;
	cout << "period:" << T << endl;

	return 0;
}

当然还有直接使用dsu的解法：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
using ll = long long;

ll A, C, M;

inline ll f(ll x)
{
    return (A * x + C) % M;
}

signed main()
{
    cin >> A >> C >> M;
    vector<ll> dsu(M, -1);
    ll x=f(0), s=0;
    while(1)
    {
        if(dsu[x]!=-1)
        {
            ll t=dsu[x];     // 预周期
            ll T=s-dsu[x];   // 周期

            cout << t << " " << T << endl;
            break;
        }

        dsu[x]=s;
        x=f(x);
        s++;
    }

    return 0;
}

这里解释下这个简洁的dsu为什么可以完成：

我们假设有这样一个序列，

1
2
3
4
                            t                             s
x[0] → x[1] → x[2] → … → x[t-1] → x[t] → x[t+1] → … → x[t+T-1]
                              ↑__________________________|
							             T = s-t

这里在x[t-1]进入环，那么dsu从x[0]开始计数，数到x[t+T-1]时计数步数为s，由x=f(x)进入x[t-1]这个已经被dsu标记过，那么此时dsu[t-1]的计数步数就是预周期t，s-t就是周期。

质因数分解

对于质因数分解，我们最直接就是使用试除法，枚举不大于 $\sqrt{m}$ 的质因子，然后不断除，直到不能整除为止，但是时间复杂度是 $O(\sqrt{m})$ ，对于大数来说显然是不够的。

1
2
3
4
5
6
7
8
9
10
11
vector<int> p,ep;
void f(int x)
{
	for(int p=2; p*p<=n; p++)if(n%p==0)
	{
		int c=0;
		while(n%p==0) n/=p,c++;
		p.push_back(p), ep.push_back(c);
	}
	if(n>1) p.push_back(n), ep.push_back(1);
}

数论解法

我们考虑一个奇合数 $n$ 其能被分解成 $n=a\times b$ ，那么其必然是两个平方数的差：

这个方法我们可以先让 $x:=\sqrt n $ ，不断向上枚举x，直到 $x^2-n$ 是个完全平方数为止，那么 $y=\sqrt{x^2-n}$ ，然后 $a=x+y$ ， $b=x-y$ ，就可以得到 $n$ 的一个因数分解。
然后继续按照上面的逻辑分解a和b，直到分解到质数为止。

这种方法对于n的两个因子十分接近时，效率会很高。

但是反之，如果n的两个因子差距很大，我们需要枚举到很大的x才能找到这个完全平方数，我们可以考虑同余式：

那么我们只要找到 $n | x^2-y^2$ 即可，或者构造 $x^2 \equiv y^2 \pmod n$ 的解，比如随机或者构造平方。

我们考虑利用最大公约数提取因子，如果n整除 (x-y)(x+y) ，这并不代表n一定整除这其中某一个因子。而我们为了质因数分解，肯定不希望(x-y)或者(x+y)不是n的平凡因子（即因子不是1或者n本身），我们需要保证n的因子都分步在这两项中。

我们假设 $n = p \cdot q$（$p, q$ 为互异的质数），那么由CRT有：

因为 $p$ 是质数，那么 $\Z_p$ 是个域，二次方程 $x^2 - y^2 = 0$ 只有两个解（域上的多项式性质）：

同理，对于模 $q$ 也有：

也就是在模n的意义下这四个解：

  • 平凡解
  • 非平凡解

因此我们对于 $n | x^2-y^2$ 再加以非平凡解限定的条件 $x \not \equiv \pm y \pmod n$ ，此时n的一部分因子必然在左右两项都有分布。

具体的分解过程如下：

  • 寻找，找到一对 $(x, y)$ 使得 $x^2 \equiv y^2 \pmod n$
  • 提取，计算 $g = \gcd(x-y, n)$
  • 拆分， $n := g \times (\frac{n}{g})$。我们将 $a$ 和 $b = \frac{n}{a}$ 视为新的待分解数
  • 终止，如果某个数通过 Miller-Rabin 等算法检测为质数，则停止对该分支的迭代

即产生了一颗递归树——因子分解树。

程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
using ll = long long;

ll mul(ll a, ll b, ll mod)
{
    return (__int128)a * b % mod;
}

ll qpow(ll a, ll b, ll mod)
{
    ll res=1;
    while(b)
    {
        if(b&1) res=mul(res, a, mod);
        a=mul(a, a, mod);
        b >>= 1;
    }
    return res;
}

bool MillerRabin(ll n)
{
    if (n < 2) return false;
    for (ll p : {2, 3, 5, 7, 11, 13, 17, 19, 23})
    {
        if (n % p == 0) return n == p;
    }

    ll d = n - 1, s = 0;
    while ((d & 1) == 0)
    {
        d >>= 1;
        s++;
    }

    auto check = [&](ll a)
    {
        ll x = qpow(a, d, n);
        if (x == 1 || x == n - 1) return true;
        for (int i = 1; i < s; i++)
        {
            x = mul(x, x, n);
            if (x == n - 1) return true;
        }
        return false;
    };

    for (ll a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022})
    {
        if (a % n == 0) continue;
        if (!check(a)) return false;
    }

    return true;
}

ll PollardRho(ll n)
{
    if(~n&1) return 2;
    while(1)
    {
        ll x=rand()%(n-2)+2;
        ll y=x;
        ll c=rand()%(n-1)+1;

        auto f=[&](ll x)
        {
            return (mul(x, x, n)+c)%n;
        };

        ll d=1;

        while(d==1)
        {
            x=f(x);
            y=f(f(y));
            d=__gcd(abs(x - y), n);
        }
        if(d!=n) return d;
    }
}

map<ll,int> mp; // <p,ep>

// 质因数分解
void factor(ll n)
{
    if(n==1) return;

    if(MillerRabin(n))
    {
        mp[n]++;
        return;
    }

	// 用PollardRho分解
    ll d=PollardRho(n);
    factor(d);
    factor(n/d);
}

图论解法