【前缀和 差分数组 数论】P6042 「ACOI2020」学园祭|省选-

本文涉及的基础知识点

C++算法：前缀和、前缀乘积、前缀异或的原理、源码及测试用例 包括课程视频
C++差分数组
数论：质数、最大公约数、菲蜀定理

「ACOI2020」学园祭

题目背景

秋天，是学习之秋，食欲之秋，更是，学园祭之秋！随着时间流逝，学园祭也越来越近。终于等到这一天，可是没想到在冲绳岛上邂逅到女装的渚同学的勇次竟然来了！中村 莉櫻（Nakamura Rio）见到这个情况，忙给渚同学换上女装。没办法，勇次已经来了，于是渚同学鼓起勇气迈出了第一步。（为什么自顾自地加提示框啊喂！）

题目描述

莉櫻为了利用这个人傻钱多的少爷，尽全力提高消费额，努力地暗示渚同学。没办法，于是渚同学想了一下，提出了一个问题：

给出一个 $n$，定义：
$$\Gamma (0)=1,\Gamma (n)=n!$$

$$A_i^j=\frac{\Gamma (i)}{\Gamma (j)}$$
求
$$\sum _{i=1}^n\sum _{j=1}^i\sum _{k=1}^j\gcd (A_{i-j}^j\times \Gamma (j),A_{j-k}^k\times \Gamma (k))$$

渚同学念着莉櫻举起的对话板上写的字：如果不能在规定时间回答出问题的话，就要把菜单全部买一遍哦！

尽管勇次钱多，但是他并不想吃得太多，因为这个问题有 $T$ 个小问题！

由于答案可能太大，请将答案对 $10086001$ 取模。

输入格式

本题有多组数据。

第一行一个整数 $T$，表示数据组数。

对于每组数据：

只有一行一个整数 $n$。

输出格式

对于每组数据，一行一个整数，表示问题的答案对 $10086001$ 取模后的值。

样例 #1

样例输入 #1

5
1
2
3
4
5

样例输出 #1

1
4
10
20
36

提示

数据范围

本题采用捆绑测试。

• Subtask 1（20 points）：$T\le 10^3$，$n\le 10^2$。
• Subtask 2（30 points）：$T\le 10^6$，$n\le 5\times 10^3$。
• Subtask 3（50 points）：$T\le 10^6$，$n\le 10^6$。

对于 $100\%$ 的数据，$1\le T,n\le 10^6$。

差分数组 前缀和 数论

差分数组 O(nn)

对于$\forall$ i,j,k，可以化简为gcd((i-j)!,(j-k)!)。
对于$\forall$j, i-j$\in$[0,n-j]，j-k $\in [0,j-1]对于$\forall$j,
gcd为0!的数量为：(n-j)+(j-1)+1=n,参数1为0，参数2取值[1,j-1]；参数2为0，参数1取值[1,n-j];参数1和参数皆取0。
gcd为1!的数量为：n-2
$⋮$
a[i]记录gcd为i!的数量 = n-2x x <= min(n-j,j-1)

枚举j:
情况一：如果 n-j <= j-1 $⟺$ n+1 <= 2j
如果n+1是偶数： $iff$ j >= （n+1）/2
如果n+1是奇数：$⟺$ j >= (n+1+1)/2
两者可以统一为： j >= j1= (n+2)/2
vDiff[0] += 1;
vDiff[n-j+1] -= 1
如果j < j1
vDiff[0] += 1;
vDiff[j] -= 1
c是vDiff的前缀和
ans +=x! *c[x]*a[x]
时间复杂度：O(n)

前缀和

如果提前预处理所有n，时间复杂度O(nn)。如果每次查询都处理，时间复杂度O(Tn)，都过不了第三个子任务。

前缀和

e[x] 记录 j < j1对c[x]的贡献
x <= j-1，j-1取值范围[0,j1-2]
e[x] = max(0,j1-2-x+1)=max(0,j1-x-1)
f[x]记录 j >= j1对c[x]的贡献
x<=n-j ,n-j的取值范围是[0,n-j1]
f[x] = max(0,n-j1-x+1)
c[x] = e[x]+f[x]
时间复杂度和差分数组相同。

数学

当n较大时,e[x]和f[x]都大于0，则c[x] = j1-x-1+n-j1-x+1 = n - 2x
a[x]c[x] = (n-2x)(n-2x)
n至少多大，e[x]和f[x]都大于0
e[x] >0 $⟺$ (n+2)/2-x-1 >= 0 $⟺$ n/2 >=x
f[x] > 0 $⟺$ n - (n+2)/2-x+1 >= 0 $iff$ n-n/2>=x ,n无论是奇数，还是偶数，都$⟺$ x <=n/2。
总结
$$c[x]=\{\begin{array}{llc}(n-2x)\ast (n-2x）& & x<=n/2\\ 0& & other\end{array}$$
g(n)记录本题答案,则 g(0)=0，g(1)=1
如果n是奇数，g(n-1)和g(n)的项数相同；如果n是偶数,g(n)比g(n-1)多一项c[n/2]。
令b = (n-1)/2$\forall$ x <= b，对g(n)的贡献-g(n-1)的贡献为：
((n-2x)(n-2x)-(n-1-2x)(n-1-2x))*x!
令a = n-2x ，上式$⟺$ aa-(a-1)(a-1)=aa-aa+2a-1=2a-1
即(2n-4x-1)x!
令y +z = g(n)-g(n-1),如果n是奇数z为0，为偶数z = c[n/2],也为0。 则 y = ${\sum }_{x:0}^b(2n-4x-1)x!$
即(2n-1)$\sum$x! - 4${\sum }_{x:0}^b(xx!)$

通过g(n-1)计算g(n)的时间负重度为：O(1)
故预处理所有的n，时间复杂度为：n。

代码

1e6个整数和换行，用cout就超时,估计用时300ms，换成printf就正常。

O(nn)

class Solution_2 {typedef  C1097Int<10086001> BI;public:vector<int> Ans(vector<int>& a) {static auto vInit = Init();vector<int> ans;for (const auto& i : a) {ans.emplace_back(vInit[i]);}return ans;}vector<int> Init() {const int N = 5'000;vector<BI> fac(N + 1, 1);vector<int> ret = { 0 };for (int i = 1; i <= N; i++) {fac[i] = fac[i - 1] * i;					BI b1 = 0;for (int x = 0; x <= i / 2; x++) {b1 += fac[x]*(i - x * 2) * (i - x * 2);}ret.emplace_back(b1.ToInt());}return ret;}};

代码

#include <iostream>
#include <sstream>
#include <vector>
#include<map>
#include<unordered_map>
#include<set>
#include<unordered_set>
#include<string>
#include<algorithm>
#include<functional>
#include<queue>
#include <stack>
#include<iomanip>
#include<numeric>
#include <math.h>
#include <climits>
#include<assert.h>
#include<cstring>#include <bitset>
using namespace std;template<class T1, class T2>
std::istream& operator >> (std::istream& in, pair<T1, T2>& pr) {in >> pr.first >> pr.second;return in;
}template<class T1, class T2, class T3 >
std::istream& operator >> (std::istream& in, tuple<T1, T2, T3>& t) {in >> get<0>(t) >> get<1>(t) >> get<2>(t) ;return in;
}template<class T1, class T2, class T3, class T4 >
std::istream& operator >> (std::istream& in, tuple<T1, T2, T3, T4>& t) {in >> get<0>(t) >> get<1>(t) >> get<2>(t) >> get<3>(t);return in;
}template<class T = int>
vector<T> Read() {int n;scanf("%d", &n);vector<T> ret(n);for(int i=0;i < n ;i++) {cin >> ret[i];}return ret;
}template<class T = int>
vector<T> Read(int n) {vector<T> ret(n);for (int i = 0; i < n; i++) {cin >> ret[i];}return ret;
}template<int MOD = 1000000007>
class C1097Int
{
public:C1097Int(long long llData = 0) :m_iData(llData% MOD){}C1097Int  operator+(const C1097Int& o)const{return C1097Int(((long long)m_iData + o.m_iData) % MOD);}C1097Int& operator+=(const C1097Int& o){m_iData = ((long long)m_iData + o.m_iData) % MOD;return *this;}C1097Int& operator-=(const C1097Int& o){m_iData = (m_iData + MOD - o.m_iData) % MOD;return *this;}C1097Int  operator-(const C1097Int& o){return C1097Int((m_iData + MOD - o.m_iData) % MOD);}C1097Int  operator*(const C1097Int& o)const{return((long long)m_iData * o.m_iData) % MOD;}C1097Int& operator*=(const C1097Int& o){m_iData = ((long long)m_iData * o.m_iData) % MOD;return *this;}C1097Int  operator/(const C1097Int& o)const{return *this * o.PowNegative1();}C1097Int& operator/=(const C1097Int& o){*this /= o.PowNegative1();return *this;}bool operator==(const C1097Int& o)const{return m_iData == o.m_iData;}bool operator<(const C1097Int& o)const{return m_iData < o.m_iData;}C1097Int pow(long long n)const{C1097Int iRet = 1, iCur = *this;while (n){if (n & 1){iRet *= iCur;}iCur *= iCur;n >>= 1;}return iRet;}C1097Int PowNegative1()const{return pow(MOD - 2);}int ToInt()const{return (m_iData + MOD) % MOD;}
private:int m_iData = 0;;
};class Solution {typedef  C1097Int<10086001> BI;public:vector<int> Ans(vector<int>& a) {static auto vInit = Init();vector<int> ans;for (const auto& i : a) {ans.emplace_back(vInit[i]);}return ans;}vector<int> Init() {const int N = 1000'000;vector<BI> fac(N + 1, 1);vector<BI> v(N + 1);vector<BI>s1(N / 2 + 1,1), s2(N / 2 + 1,0);for (int x = 1; x <= N / 2; x++) {fac[x] = fac[x - 1] * x;s1[x] = s1[x - 1] + fac[x];s2[x] = s2[x - 1] + fac[x] * x;}vector<int> ret = { 0 };for (int n = 1; n <= N; n++) {	const int b = (n - 1) / 2;v[n] = v[n - 1]+ s1[b]*(2*n-1) -s2[b]*4;ret.emplace_back(v[n].ToInt());}return ret;}};int main() {
#ifdef _DEBUGfreopen("a.in", "r", stdin);
#endif // DEBUG		auto a = Read<int>();	
#ifdef _DEBUG	//printf("N=%d", n);Out(a, ",a=");//Out(b, ",b=");
#endif	auto res = Solution().Ans(a);for (const auto& i : res) {printf("%d\r\n",i);}return 0;
}

单元测试

	vector<int> a;TEST_METHOD(TestMethod11){a = { 1,2,3,4,5 };auto res = Solution().Ans(a);AssertEx({ 1,4,10,20,36 }, res);}TEST_METHOD(TestMethod12){a.clear();for (int i = 1; i <= 5000; i++) {a.emplace_back(i);}auto res1 = Solution_2().Ans(a);auto res = Solution().Ans(a);AssertEx(res1, res);}

扩展阅读

我想对大家说的话
工作中遇到的问题，可以按类别查阅鄙人的算法文章，请点击《算法与数据汇总》。
学习算法：按章节学习《喜缺全书算法册》，大量的题目和测试用例，打包下载。重视操作
有效学习：明确的目标 及时的反馈 拉伸区（难度合适） 专注
闻缺陷则喜(喜缺)是一个美好的愿望，早发现问题，早修改问题，给老板节约钱。
子墨子言之：事无终始，无务多业。也就是我们常说的专业的人做专业的事。
如果程序是一条龙，那算法就是他的是睛
失败+反思=成功 成功+反思=成功

视频课程

先学简单的课程，请移步CSDN学院，听白银讲师（也就是鄙人）的讲解。
https://edu.csdn.net/course/detail/38771
如何你想快速形成战斗了，为老板分忧，请学习C#入职培训、C++入职培训等课程
https://edu.csdn.net/lecturer/6176

测试环境

操作系统：win7 开发环境： VS2019 C++17
或者 操作系统：win10 开发环境： VS2022 C++17
如无特殊说明，本算法用**C++**实现。