jzoj6748 回文串（字符串结论计数，反演, min25）

题面

求长度不超过n的双回文串个数。
$n\le 10^9$

分析

这是个结论题，主要难点是分析出弱回文串的结构。（或者打表发现他）

弱双回文串=第一个回文串可为空的双回文串。
最小弱回文指的是，最小整周期为n的长度为n的弱回文串。

1. 一个有多个划分方法的弱双回文串=某个最小弱双回文串的若干复制
2. 弱双回文串的划分方案数即为，最大的可能复制次数。
3. 双回文串个数=弱双回文串-最小回文串

对拍，上述结论均正确，开始你的数论表演：
设：
$f(n)$为长度为n的最小弱双回文串个数，
$g(n)$为长度为n的弱双回文划分个数，
$h(n)$为长度为n的弱双回文个数。

$g(n)$只需分类讨论一下奇偶，可得
$g(n)=\{\begin{array}{l}n{C}^{(n+1)/2},n\%2=1\\ n/2(C^{n/2}+C^{n/2+1}),n\%2=0\end{array}$

由上述，$$g(n)=\sum _{d|n}f(d)\times (n/d)$$，反演得$$f(n)=\sum _{d|n}g(d)(n/d)\mu (n/d)$$
并且有$$h(n)=\sum _{d|n}f(d)$$
接下来是求解${\sum }_{i=1}^{n}h(n)$
$$=\sum _{i=1}^n\sum _{d|i}f(d)$$
$$=\sum _{d=1}^{n}f(d)(n/d)$$
$$=\sum _{d=1}^n(n/d)\sum _{k|d}g(k)(d/k)\mu (d/k)$$
$$=\sum _{k=1}^{n}g(k)\sum _{d^{\prime }=1}^{n/k}(\frac{n/k}{d^{\prime }})\mu (d^{\prime })d^{\prime }$$
该式可以分块，$g$是个等差比可以快速求和。现在要快速算后面的关于n/k的函数。
$$T(n)=\sum _{i=1}^n(\frac{n}{i})\mu (i)i$$，注意括号里是下取整。
$$=\sum _{j=1}^n\sum _{d|j}\mu (d)d$$，这是个积性函数前缀和，直接min25即可。（也可以杜教筛）

接下来计算最小回文数，基本同理：
设$b(n)$为长度为n的回文数，$c(n)$为长度为n的最小回文。注意到$b(n)$很好算。
$$b(n)=\sum _{d|n}c(d)$$
反演，$$c(n)=\sum _{d|n}b(d)\mu (n/d)$$
对$c(n)$求前缀和$$\sum _{i=1}^{n}c(i)$$
$$=\sum _{i=1}^n\sum _{d|i}b(d)\mu (i/d)$$
$$=\sum _{d=1}^{n}b(d)\sum _{i=1}^{n/d}\mu (i)$$
分块+min25/杜教筛即可。

结论的证明

（由于我的脑抽，以下弱回文=弱双回文
接下来是对结论的证明，学艺不精可能有伪证。。。请各位抱着批判的眼光看待。。。
首先有弱周期引理：若有周期$p,q$且$p+q\le |S|$，则$gcd(p,q)$为S的一个周期。

Lemma 2

一个有两个以上划分的弱回文是整周期串。
只需找出两个周期$a,b$使得$gcd(a,b)|s$。设两个划分的第一个回文长度差值为$t$，容易发现s同时有$t$的周期和border.
t的周期：大回文前缀包含小回文前缀，会产生长度为他们的差的周期。
t的border：通过一些回文翻转即可发现。

Lemma 3

弱回文串$S$的最小整周期是最小弱回文串，并且该串的弱回文划分数恰好为$|S|/最小整周期$。

设某个弱回文串的最小整周期所对应的串为$T$，由$Lemma2$，T至多只有一个弱回文划分。
下面证明$T$存在一个弱回文划分。
当$T=S$时引理成立，否则设串$S=pq$为其弱回文划分，取包含至少一个T的部分，不失一般性地我们取$q$，$q$是若干个T与一个T的可空前缀拼接而成。由于q是回文串，容易发现该前缀是回文串。并且去掉该前缀之后剩下的$T^{\prime }$也是回文串。这就证明了$T$是一个最小弱回文串（只存在唯一划分）。

然后考虑其弱回文划分数。首先至少存在$|S|/|T|$个。然后，由上述分析我们回文划分的分界点必须取在$T$的回文划分分界点上，否则会与$T$只存在唯一划分矛盾。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mo = 1e9 + 7, N = 1e6 + 10;

int n, c;
int is[N], mu[N];
ll p[N];
int loc[N], loc2[N], w[N], B, precnt[N], presum[N], preT[N];
#define id(x) ((x) <= B ? loc[x] : loc2[n / (x)])
void init(int n) {
	mu[1] = 1;
	for(int i = 2; i <= n; i++) {
		if (!is[i]) p[++*p] = i, mu[i] = -1;
		for(int j = 1; p[j] * i <= n; j++) {
			is[p[j] * i] = 1;
			if (i % p[j] == 0) {
				break;
			} else {
				mu[i * p[j]] = -mu[i];
			}
		}
		precnt[i] = precnt[i - 1] + (is[i] == 0 ? 1 : 0);
		presum[i] = (presum[i - 1] + (is[i] == 0 ? i : 0)) % mo;
		preT[i] = (preT[i - 1] + (is[i] == 0 ? 1 - i : 0)) % mo;
	}
}

ll cnt[N], sum[N], smu[N], aT[N];



void min25() {
	B = sqrt(n);
	for(int i = 1; i <= n; ) {
		int v = n / i, r = n / v;
		w[++(*w)] = n / i;
		if (v <= B) loc[v] = *w; else loc2[i] = *w;
		cnt[*w] = v - 1;
		sum[*w] = (2ll + v) * (v - 1) / 2 % mo;
		i = r + 1;
	}
	for(int i = 1; p[i] * p[i] <= n; i++) {
		for(int j = 1; w[j] >= p[i] * p[i]; j++) {
			cnt[j] -= cnt[id(w[j] / p[i])] - precnt[p[i] - 1];
			sum[j] = (sum[j] - p[i] * (sum[id(w[j] / p[i])] - presum[p[i] - 1])) % mo;
		}
	}
	
	for(int i = 1; i <= *w; i++) {
		smu[i] = -cnt[i];
		aT[i] = cnt[i] - sum[i];
	}
	
	for(int i = p[0]; i; i--) {
		for(int j = 1; w[j] >= p[i] * p[i]; j++) {
			for(int e = 1, b = p[i]; b * p[i] <= w[j]; e++, b *= p[i]) {
				if (e == 1)
					smu[j] = (smu[j] + (-1) * (smu[id(w[j] / b)] - (-1) * precnt[p[i]]));
				aT[j] = (aT[j] + (1 - p[i]) * (aT[id(w[j] / b)] - preT[p[i]]) + (1 - p[i])) % mo;
			}
		}
	}
}

ll ans;

ll ksm(ll x, ll y) {
	ll ret = 1; for(; y; y >>= 1) {
		if (y & 1) ret = ret * x % mo;
		x = x * x % mo;
	}
	return ret;
}

ll g(ll n) {
	if (n & 1) return n * ksm(c, (n + 1) / 2);
	else return n / 2 * (ksm(c, n >> 1) + ksm(c, (n >> 1) + 1)) % mo;
}

ll A1, B1, A2, B2;
ll _sumg(ll n) {
	ll p1 = ksm(c, n / 2 + 1) * (A1 * (n / 2) % mo + B1) - c * B1;//odd
	p1 %= mo;
	ll dt = 0;
	
	ll p2 = ksm(c, n / 2 + 1) * (A2 * (n / 2) % mo + B2) - c * B2; //even1
	p2 %= mo;
	
	p2 = p2 * (1 + c) % mo; //even2
	return (p1 + p2) % mo;
}

ll sumg(ll n) {
	ll rv = _sumg(n - (n & 1)) + ((n & 1) ? g(n) : 0);
	return rv % mo;
}

ll b(ll n) {
	return ksm(c, (n + 1) >> 1);
}

ll _sumb(ll n) {
	ll cnt = n >> 1;
	return 2 * (ksm(c, cnt + 1) - c) % mo * ksm(c - 1, mo - 2) % mo;
}

ll sumb(ll n) {
	//ll ret = 0;
	ll rv = _sumb(n - (n & 1)) + (n & 1 ? b(n) : 0);
	//for(int i = 1; i <= n; i++) ret = (ret + b(i)) % mo;
	//assert((ret + mo) % mo == (rv + mo) % mo);
	return rv;
}

ll sumu(ll x) {
	return smu[id(x)] + 1;
}

ll T(ll x) {
	return aT[id(x)] + 1;
}

int main() {
	freopen("string.in","r",stdin);
	freopen("string.out","w",stdout);
	cin >> n >> c; init(1e6); min25();
	A1 = 2 * ksm(c - 1, mo - 2) % mo;
	B1 = - (1 + A1) * ksm(c - 1, mo - 2) % mo;
	A2 = 1 * ksm(c - 1, mo - 2);
	B2 = - A2 * ksm(c - 1, mo - 2) % mo;
	for(int l = 1; l <= n; ) {
		int r = n / (n / l);
		ans = (ans + (sumg(r) - sumg(l - 1)) * T(n / l)) % mo;
		ans = (ans - (sumb(r) - sumb(l - 1)) * sumu(n / l)) % mo;
		l = r + 1;
	}
	cout << (ans + mo) % mo << endl;
}