数学问题总结
数学问题总结
构造
-
题目大意：平面直角坐标系里面有很多点，都是整数点，求一条直线，使得在直线上以及上方的点与在直线下方的点的个数比为\(a:b\)，其中点的坐标值的绝对值均\(\leq10^5\)，所求直线的斜率应满足\(1\leq k\leq10^{12}\) 。题目保证\(ab\neq0\) 且\((a+b)|n\) 。
在做这道题的时候，我发现了坐标值有限，并且斜率可以很大，所以考虑了斜率很大的情况，当斜率很大时，差不多就是一条竖线，则大概可以只关心横坐标的情况，所以在点斜式中我直接令\(k=10^{12},y_0=0\)，然后二分求\(x_0\)。但是这样存在问题，一个是我们得关心点的纵坐标是大于等于0还是小于0，前者应算在上面，后者是下面，第二个是取\(y_0=0\) 不一定有解。于是就有了题解的解法。
我们把点按照横坐标进行排序，当横坐标相同时，按照纵坐标降序排序，这样的话，考虑令直线斜率很大，且经过第\(\frac{na}{a+b}\)个点，则我们就构造出来了答案。
公式化简
-
题目大意：见原题。
这个题充分考察了数学式子的化简，以及前缀和的使用，和某逆序对问题的解题关键类似，都是式子化简。本题我还有一些化简感觉不是很自然，先占坑。
能够直接把最终答案的式子写出来的题，一般都要进行公式化简，因为朴素求值的话复杂度会不对，要充分利用之前计算的结果进行优化，比如前缀和。
如本题，\(ans=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{[\frac{n-i}{2}]}(2i+2j)(num[i]+num[i+2j])(color[i]==color[i+2j])\) ，朴素计算复杂度\(O(n^2)\) 。
我们发现这个式子有很烦人的地方，比如那个条件式，因为有这个东西，导致我们一点点东西都没法从内层循环中提取出来。如果计算式中所有的都默认满足这个条件就好了。为了保证这一点，我们需要考虑把满足这个条件的格子都放在一堆中，这样的话整个纸带上的格子就被分成了若干个堆，其中每个堆里的格子两两满足这个条件。这个条件怎么筛出来呢？它蕴含两点：一个是颜色相同，一个是下标满足那个关系，可以看做等差数列，也可以看做奇偶性相同。这样的话，分堆问题容易解决了：按照颜色和奇偶性两个指标即可分类。考虑某个堆，显然堆里面任意两个元素组合都可以做出贡献，且所有的贡献就是他们的和。记这个堆里面有\(cnt\)个元素，分别为\(i_1,i_2……i_{cnt}\)，则最终的贡献是\((i_1+i_2)(num[i_1]+num[i_2])+(i_1+i_3)(num[i_1]+num[i_3])+……\)\(+(i_1+i_{cnt}(num[i_1]+num[i_{cnt}]))+……\)\(+(i_{cnt}+i_{cnt-1})(num[i_{cnt}]+num[i_{cnt-1}])\) 。显然朴素统计是\(O(cnt^2)\)的复杂度。考虑化简，合并同类项与构造前缀和：先只看\(i_1\)相关的那部分，最后整理的结果应该是\(i_1(num[i_2]+num[i_3]+…num[i_{cnt}])+(cnt-1)(i_1num[i_1])+(i_2+i_3+…\)\(+i_{cnt})num[i_1]+(i_2num[i_2]+i_3num[i_3]+…+i_{cnt}num[i_{cnt}])\) 。显然括号中涉及到的和式都可以优化，用前缀和就能搞定。
位运算相关
-
题目大意，有n本书，这n本书其实是好几种书，其中保证每种书的本书是k的倍数，现在我往里面再插入一本书，请求出来我插的那本书是哪种书。
我们知道，k为偶数的时候，就是一道经典的面试题，用到的是异或的思想，即由于插入了一本书，所以肯定有某种书为奇数本了。而两本相同种类的书，异或之后结果是0，所以，把所有书的种类做一个异或，最后得到的数就是新插入的书的种类。
那么k为奇数的时候呢？这里要介绍一下异或运算的本质：不进位加法。0^0=0 0^1=1 1^0=1 1^1=0，是不是？裸的异或是二进制运算，也就是说，两个一样的做异或运算，得到的是0，那么我们可以思考，能不能发明一种k进制异或运算，使得有k个一样的数进行这个运算时，结果是0呢？如果理解了不进位加法的思想，会发现这是可行的。我们把书的种类变成k进制数，然后把k的倍数本书的种类的k进制表示都做不进位加法（最后统一对k进制位取模）,那么，如果这些书编号都是一样的，那么他们的k进制表示也是一样的，所以它们的每个k进制位加起来之后肯定都是k的倍数（因为是k的倍数本一样的书做不进位加法）,每一位都对k取模之后肯定是0啦。所以，我们对所有的n+1本书都做k进制不进位加法，然后取模，最后把得到k进制数变回10进制，就是插入的书的种类啦。至于复杂度，如果种类最大是A，则复杂度为\(O(nlog(A))\)。当然，这个题这样还过不了，因为出题人太毒瘤了，把A设得特别大，所以这道题目正解是把64位数变成4个16位数（2进制），然后对于n+1本书，看看每本书的这4个部分出现了多少次，把不是k的倍数的部分拎出来，最后一整合就是插入的书了。
不进位加法思路代码如下(75pts)：
#include <bits/stdc++.h> #define ll long long using namespace std; const int N=1e7+8; ll n,k,ans; ll tmp[66]; int main(){ ll x; scanf("%lld %lld",&n,&k); if(k%2==0){ scanf("%lld",&ans); for(int i=2;i<=n;i++){ scanf("%lld",&x); ans=ans^x; } printf("%lld\n",ans); } else { for(int i=1;i<=n;i++){ scanf("%lld",&x); for(int j=0;x;j++,x=x/k) tmp[j]+=x%k; } for(int j=0;j<64;j++) tmp[j]%=k; for(int j=63;j>=0;j--) ans=ans*k+tmp[j]; printf("%lld\n",ans); } return 0; }正解代码如下(不需要读入优化)：
#include <bits/stdc++.h> #define ll long long #define gc pa==pb&&(pb=(pa=buf)+fread(buf,1,100000,stdin),pa==pb)?EOF:*pa++ using namespace std; ll n,k,ans; ll tmp[5][(1<<16)+9]; static char buf[100000],*pa(buf),*pb(buf); inline ll readint(); int main(){ ll x; n=readint(); k=readint(); // scanf("%lld %lld",&n,&k); for(int i=0;i<n;i++){ x=readint(); // scanf("%lld",&x); tmp[3][(x>>48)&0xffff]++; //统计每个数的块出现的个数 tmp[2][(x>>32)&0xffff]++; //统计每个数的块出现的个数 tmp[1][(x>>16)&0xffff]++; //统计每个数的块出现的个数 tmp[0][x&0xffff]++; //统计每个数的块出现的个数 } for(ll i=3;i>=0;i--){ for(ll j=0;j<65536;j++){ if(tmp[i][j]%k){ ans|=(j<<(i*16)); break; } } } printf("%lld\n",ans); return 0; } inline ll readint() { ll x=0;char c=gc; while(c<'0'||c>'9')c=gc; for(;c>='0'&&c<='9';c=gc)x=x*10+(c&15); return x; }
数论
-
题目大意：从数1-n个中取出k个数，使得这k个数最大公因数最大，求这个最大公因数。
可能比较容易想到的一个方法是筛法：我可以尝试从大到小枚举这个最大公因数，然后看看能筛去哪些数，第一次枚举到的能够筛去超过k个数的数就是最大公因数了。不过这个题n的规模是1e9，会TLE。这时候我们要换方法，毕竟数论题有些就是结论嘛，说不定有一些结论可以用到。其实有了刚才筛法的经验，或许我们不难想到，要想用筛法筛出来k个数，那么这个最大公约数肯定不能大于\([\frac{n}{k}]\)，不然根本筛不够k个，即\(gcd<=[\frac{n}{k}]\)。另一方面，如果我们求出了\(gcd\)，那么肯定是筛去了\(gcd,2gcd...kgcd\) ，如果\(gcd=[\frac{n}{k}]\)，那么显然有\(kgcd<=n\)，即\(gcd=[\frac{n}{k}]\) 是合法的。所以，要求的最大的\(gcd\)就是\(\frac{n}{k}\)了。虽然我们不一定能直接想到这个结论，但是，我们还是能够通过筛法去发现这个结论的。
-
题目大意：对于代数式\(ax+by\)，已知\(a,b\) 为正整数，由输入给定，且令\(x,y\)取遍所有非负数值，并保证这时\(ax+by\)不能组合出来的整数\(c\)有最大值，求这个最大的\(c\)
朴素的想法是，我们可以发现\(c\)的最小可能取值是\(min(a,b)-1\) ，最大应该不会超过\(ab\) （显然，取\(x=b,y=0\)能使得代数式组合出\(ab\)，所以肯定不超过\(ab\) ），所以我们从大到小枚举\(c\)（答案没有单调性，只能顺序枚举），每枚举一次，就做一次\(exgcd\)来判断一下，先看有没有整数解，再看有没有正整数解，第一个不合理的\(c\) 就是答案了。这个思路可以轻易通过\(60\%\) 的数据。
剩下的\(40\%\)的数据，\(a,b\)很大，原算法会超时。由于我们前面的算法是基于枚举的，所以可能会比较慢，但\(exgcd\)这个工具应该还是有必要用的，只不过需要变形一下。根据题意，如果\(c\)求出来了，那么令\(d=c+1\)，\(ax+by=d\)是有合理解的。那么，我们就是想求这样一个最大的\(d\)：\(ax+by=d\)有合理解，且\(ax+by=d-1\)没有合理解(这一步非常重要，是数学中的一种很常见的思路，化归的思想)。如果\(ax+by=d-1\)也有合理解，至少有整数解，假设\(ax_1+by_1=d\)，那么就有\(a(x_1-\Delta x)+b(y_1-\Delta y)=d-1\) ,化简一下并带入上式有\(a\Delta x+b\Delta y=1\) 。要解出来整数\(\Delta x\)和\(\Delta y\)，需要\(gcd(a,b)=1\)。也就是说如果没有\(gcd(a,b)=1\)的话，这个题目的数据就有锅，就没法保证有最大的\(c\)。假设题目数据没锅，那么能解出来可能的\(\Delta x\)和\(\Delta y\) 。然后就不会分析了。
偷偷看了一眼提示，由于\(gcd(a,b)=1\)，所以肯定有整数解，我们现在求出两组解，一组是\((x',y')\)，另一组是\((x'',y'')\)，其中，\(x'\)是这个方程的整数解中\(\Delta x\)非负且最小的，\(y''\)是这个方程整数解中\(\Delta y\)非负且最小的。 要想让\(ax+by=d-1\)没有合理解，需要\(x_1-\Delta x\)和\(y_1-\Delta y\)中至少有一个恒为负数。如果\(\Delta x\)取\(x'\),都能让\(x_1-\Delta x=x_1-x'<0\) ，或者\(\Delta y\)取\(y''\)，都能让\(y_1-\Delta y=y_1-y''<0\) ，二者中至少有一个成立，那就肯定不会有非负整数解了。为了让\(d\)尽可能小，我们取\(x_1=x'-1\)，\(y_1=y''-1\) ，则\(a(x'-1)+b(y''-1)\)就是我们所说的最大的\(d\)，\(d-1\)就是最大的不能被组合出来的数。就这样，答案就被构造出来了。虽然这个比某个很暴力一个式子就能艹过去的方法要弱，但是，能摆脱大家对于小学奥数的阴影。
小学奥数做法：如果\(a,b\)都至少选一个的话，那么最大的\(ax+by\)拼不出来的是\(ab\),怎么来的？大胆猜想，无需证明。那么如果\(a,b\)可以不选的时候，只要在\(ab\)的基础上减去\(a+b\)即可了。严谨的证明这里我不再给出。
我认为的标准做法代码如下：
#include <bits/stdc++.h> #define ll long long using namespace std; void exgcd(ll a,ll b,ll c,ll &x,ll &y); ll lcm(ll a,ll b); ll gcd(ll a,ll b); int main(){ ll a,b,x,y,minx,miny,k,l; scanf("%lld %lld",&a,&b); k=lcm(a,b)/a; l=lcm(a,b)/b; exgcd(a,b,1,x,y); minx=((x%k)+k)%k; miny=((y%l)+l)%l; printf("%lld\n",a*(minx-1)+b*(miny-1)-1); return 0; } ll lcm(ll a,ll b){ return (a*b)/gcd(a,b); } ll gcd(ll a,ll b){ if(b==0){ return a; } return gcd(b,a%b); } void exgcd(ll a,ll b,ll c,ll &x,ll &y){ if(b==0){ x=c; y=0; return; } exgcd(b,a%b,c,x,y); ll tmp=x; x=y; y=tmp-(a/b)*y; return; } -
题目大意：找到n的第k小的因子（不一定是素因子），或者告知其不存在。n的规模是1e15，k的规模是1e9
根据数据范围可以知道，这道题起码要是\(O(\sqrt n)\)的算法，或者找规律\(O(1)\) 。
如果n 有第k小的因子的话，那么要么这个因子小于等于\(\sqrt n\)，要么大于\(\sqrt n\)。我们可以枚举从1到\(\sqrt n\)的所有的数，判断是不是因子，如果在这个范围内就找到了，那么皆大欢喜，如果没找到，假如在小于等于\(\sqrt n\)的情况下找到了m个因子，如果\(\sqrt n\)也是n的因子的话，那么一共就有\(2m-1\)个因子。如果\(2m-1<k\)，则无解，否则，我们可以根据对称性，首先我们知道第k个因子以及后面的因子一共有\(2m-1-k+1=2m-k\)个，即从右往左数，k是第\(2m-k\)个数，所以在左边我们直接找第\(2m-k\)个，然后被n除就好了。如果\(\sqrt n\)不是n的因子，那么一共有\(2m\)个因子。如果\(2m<k\)，则无解，否则根据对称性，我们要找左边的第\(2m-k+1\)个元素，被n除就是结果了。
-
题目错误大意：数列\(2,3,……n+1\) ，每个数字都要刷上一种颜色，且如果一个数是另一个数的因子，那么这两个数要刷不同的颜色，且要尽量少用颜色。给定\(n\)，求至少要用多少种颜色。
尝试了一下知道，\(n=2\)时需要1种，\(n=3……6\)时需要2种，\(n=7\)时需要3种。这样来看，似乎和取以2为底的对数有关系。考虑了一下：给定了一个\(n\)，然后我们写出来那个序列，现在我们截取序列的第\(2^i\) 到第\(2^{i+1}\)项。显然，这个区间里面不可能有倍数关系，因为即使是右端点的数，也没法大于左端点的2倍，所以更不可能是其他倍了。而如果是第\(2^i\)到第\(2^{i+1}+1\)项，那么右端点正好是左端点的2倍，就不能使得整个区间染成同一种颜色了。然后，我也不知道对不对。
然后我他妈的心态崩了！突然发现题目翻译错了！
题目正确大意：数列\(2,3,……n+1\) ，每个数字都要刷上一种颜色，且如果一个数是另一个数的质因子，那么这两个数要刷不同的颜色，且要尽量少用颜色。给定\(n\)，求至少要用多少种颜色。
纠正题意之后，因为有质因数，所以可能考虑线性筛。显然，质数全都涂一种颜色是可以的，而合数之间，由于不能有一个合数是另一个合数的质因子，所以合数也可以染相同的颜色。所以，当\(n<3\) 时，1种，\(n>=3\)时，两种。
-
题目大意：求 \(\sum\limits_{i=1}^n \gcd(i, n)\) 。
这道题n的规模是1e9，即使是线性方法都会超时，所以必须考虑根号算法或者推公式。在这里的处理方法是用欧拉函数。我们可以发现，\(\gcd(n,k)=d<=>\gcd(\frac{n}{d},\frac{k}{d})=1\) 。假如使得\(\gcd(n,k)=d\) 的\(k\) 有\(f(d)\)个，那么使得\(\gcd(\frac{n}{d},\frac{k}{d})=1\) 的\(k\)也有\(f(d)\) 个，由欧拉函数的定义知，\(f(d)=\phi(\frac{n}{d})\) 。所以可以这样考虑这个题目：枚举\(gcd\)，当\(gcd\) 取\(k\)的时候，如果\(n\% k=0\) ，则说明\(n\)有这个因数\(k\)，那么求\(f(\frac{n}{k})\)，就是求小于等于\(\frac{n}{d}\)的且和它互质的数的个数。显然，如果某个数\(i\)和\(\frac{n}{d}\)互质，那么就说明\(gcd(n,i)=d\)，所以$k*f(\frac{n}{k}) $ 就是所有 \(gcd(n,i)=k\) 的数贡献的和。顺带地，我们也可以把\(\frac{n}{k}*f(k)\)求出来，这样的话，我们只需要枚举\(gcd\)到\(\sqrt n\) 就好了。
-
题目大意：在\([1,n!]\)中，求有多少个和\(m!\)互质的数，结果对数\(r\)取模。多组数据\(n,m\)输入。
题目保证了\(m<=n\) ，所以答案组成中有一部分是\(\phi (m!)\) ，另一部分是\([m!+1,n!]\) 之间的和\(m!\)互素的数，
并且注意到\(m!\)包含了\([1,m]\)中所有的素因子。然后又注意到一件有趣的事情：当\(m>=2\)时，\(m!\)是偶数，也就是说，\([m!+1,n!]\)中所有的偶数（即一半的数）肯定都不与\(m!\)互素，剩下的奇数中，质数一定与\(m!\)互素，合数看情况，要是能快速求出一个区间之内所有和某个数互素的数就好了。
如果考虑欧拉函数单点求值公式的话，只要我们知道所有的小于等于\(m\)的素数，知道\(n!\) ，那么就可以直接求出来\(ans=n!*\Pi (1-\frac{1}{p_i})\) ，其中\(p_i\)取所有小于等于\(m\)的素数。下面考虑怎么求出来这个数。显然\(n!\)可以直接预处理+取模，\(1e7\)之内的素数可以线性筛预处理，\(p_i\)的逆元可以线性预处理， \(\Pi (1-\frac{1}{p_i})=\Pi(\frac{p_i-1}{p_i})\)在有了逆元之后可以线性预处理。都预处理好了之后，就可以\(O(1)\)询问了。
这道题给我们的启示是：如果知道一个数的素因数分解，那么就容易求出在某个区间内有多少个和它互素的数。
代码就是线性筛模板+线性求逆元模板+线性求欧拉函数模板：
#include <bits/stdc++.h> #define ll long long #define INF 9999999999 using namespace std; const int N=1e7+9; ll prime[N],inv[N],power[N],ratio[N],cnt; bool isprime[N]; ll n,m,t,mod; void getprime(); void getpower(); void getinv(); void getratio(); int main(){ scanf("%lld %lld",&t,&mod); getprime(); getpower(); getinv(); getratio(); // for(int i=1;i<50;i++){ // printf("%lld %lld %lld %lld\n",prime[i],power[i],inv[i],ratio[i]); // } while(t--){ scanf("%lld %lld",&n,&m); printf("%lld\n",(power[n]*ratio[m])%mod); } return 0; } inline void getratio(){ ratio[1]=1; for(ll i=2;i<N;i++){ if(isprime[i]){ ratio[i]=(ratio[i-1]*(((i-1)*inv[i])%mod))%mod; } else { ratio[i]=ratio[i-1]; } } } inline void getinv(){ inv[1]=1; for(ll i=2;i<N;i++){ inv[i]=((mod-mod/i)*inv[mod%i])%mod; } } inline void getpower(){ power[1]=1; for(ll i=2;i<N;i++){ power[i]=(power[i-1]*i)%mod; } } inline void getprime(){ for(ll i=2;i<N;i++){ isprime[i]=true; } for(ll i=2;i<N;i++){ if(isprime[i]){ prime[++cnt]=i; } for(ll j=1;j<=cnt && i*prime[j]<N;j++){ isprime[i*prime[j]]=false; if(i%prime[j]==0){ break; } } } } -
题目大意：见原题。
在做这道题之前，先说一下几个定理的应用。
费马小定理：若\(p\)为素数，则\(a^p\equiv a(mod\ p)\)，进一步，若\(\gcd(a,p)=1\)，则\(a^{p-1}\equiv 1(mod\ p)\) 。
欧拉定理：若\(\gcd(a,n)=1\)，则\(a^{\phi(n)}\equiv 1(mod\ n)\) 。
这两个定理都可以用来求\(a\)相对于一个模数的逆元，如果用费马小定理，则需要模数为素数；如果用欧拉定理，则需要模数与\(a\)互素。使用时，转化为快速幂计算即可。
另一个应用是，遇到求\(a^b\ mod\ n\)的时候，如果满足\(\gcd(a,n)=1\)，则可以知道\(a^{\phi(n)}\equiv 1(mod\ n)\) ，在已知这个的前提下，我们可以缩小\(b\)的范围，即让\(b=b\ mod\ \phi(n)\) ，然后就容易计算了。
本题是求\(2^{2^{2{…}}}\ mod\ p\)，指数大得不可想象，所以可能考虑用定理把指数控制一下，以求出值。每组询问会给出\(p\) ，但是不保证\(p\) 为质数。如果给出的\(p\) 能满足\(\gcd(2,p)=1\)，则可以考虑用欧拉定理，变成求\(2^{2^{2^{…}}mod\ \phi(p)}\ mod\ p\) ，这样对指数继续递归，最终\(\phi(p)\)一定会变成\(\phi(2)=1\)，这个时候，显然\(2^{2^{2^{…}}}\ mod\ 1=0\) ,在这时候就可以回溯了。但是并不是每次递归都能够让\(\gcd(2,p)=1\) 的，所以我们要考虑一个更一般的欧拉定理：
若\(b>=\phi(m)\)，则\(a^b\equiv a^{b\%φ(m)+φ(m)}(mod\ m)\) 。
并且原问题转化一下，就是求\(2^x\equiv x\ (mod\ m)\)。这道题也让我回忆了递归的易于理解的方法。
-
题目大意：求区间\([1,n]\)内所有的无序数对\((x,y)\)的个数。其中\((x,y)\)满足\(\gcd(x,y)\)为素数。\(n=1e7\)
有了前面几道题的经验，我们不难发现这道题的答案就是$\sum\limits _{i=1}^{n}\ \ \sum\limits _{k|i\ and\ k\ is \ prime}{\phi(\frac{i}{k})} $。不过可惜的是，这个没法像P2303一样利用对称性减少循环次数。
看了题解，发现此题又是另一种想法：考虑枚举素数，把题目转化为求区间\([1,\frac{n}{p}]\)内所有的\(\gcd(x',y')=1\)的无序数对的个数，其中\(x=x'p,y=y'p\)。这样的话\(O(n)\)枚举\(p\)，然后依托前缀和，\(O(1)\) 计算有多少对最大数是\(\frac{n}{p}\)的这样的数对。关于前缀和的处理，我们要线性求出欧拉函数的前缀和，需要基于线性筛，由于每个合数都是被它最小的素因子筛去的，所以可以利用欧拉函数的两个递推公式去做，即若\(\gcd(p,a)=1\)，则\(\phi(pa)=(p-1)\phi(a)\)，否则，\(\phi(pa)=p\phi(a)\)。这样就好了。
如果是有序数对呢？除了\(x=y\)的情况，我们都要加倍。所以，对于每个前缀和，我们乘2减1就好了。
代码如下，主要是依靠线性筛线性预处理欧拉函数：
#include <bits/stdc++.h> #define ll long long #define INF 999999999999 using namespace std; const int N=1e7+9; ll n,sum[N],euler[N],prime[N],cnt,ans; bool isprime[N]; void pre(); int main(){ scanf("%lld",&n); pre(); for(ll i=0;i<cnt;i++){ ans=ans+2*sum[n/prime[i]]-1; } printf("%lld\n",ans); return 0; } void pre(){ for(ll i=2;i<=n;i++){ isprime[i]=true; } euler[1]=1; for(ll i=2;i<=n;i++){ if(isprime[i]){ prime[cnt++]=i; euler[i]=i-1; } for(ll j=0;j<cnt && i*prime[j]<=n;j++){ isprime[i*prime[j]]=false; if(i%prime[j]==0){ euler[i*prime[j]]=(prime[j])*euler[i]; break; } else { euler[i*prime[j]]=(prime[j]-1)*euler[i]; } } } for(ll i=1;i<=n;i++){ sum[i]=sum[i-1]+euler[i]; } } -
题目大意：有一个数列，取出一个最长的符合要求的子列，输出子列长度以及这个子列的元素。这个子列要满足的条件是：里面所有元素的最小公倍数不超过\(m\)。
之前有一道类似的题目是：有一个集合，希望从中取出一些数，使得它们的最大公约数不是1，且取出的数尽可能多，求最多能取出来多少个数。那个题目的做法是枚举gcd，看看每次能选出来多少个数，这里是用了桶进行了优化，能够迅速知道筛掉了哪些数。
这里类似，我们枚举\([1,m]\)内的数\(y\)，看看数列中有多少数能够整除这个数。显然，对于所有的\(x|y\)，他们的最小公倍数不会超过\(y\)，更不会超过\(m\)。所以我们的方法就可以和那个题类似了。
不过这个题没有那个题友好，这样做无论是时间还是空间都是过不了的。时间上，复杂度为\(O(m\sqrt m)\) ，不行；空间上，数据范围达到了\(1e9\)，用桶开不下，不过空间问题可以用离散化来解决。
为了加速，我们考虑优化。之前是枚举的\(lcm\)，发现会超时，现在我们想到了一个类似的可以看某个数是不是某个数的约数的方法：埃氏筛。我们考虑枚举数列中的数\(x\)，然后用埃氏筛的思路，在这个数的各个倍数\(y\)上增加1的贡献，表示\(x\)是\(y\)的约数，所以选择\(y\)的时候\(x\)会有贡献。
-
题目大意：对于\(x,y\in[1,n]\)，求有多少对这样的有序对\((x,y)\)，使得\(\frac{x}{y}\)能写成有限小数。
这道题有三档部分分，\(n\)的规模一个是\(1e3\)，一个是\(1e7\)，一个是\(1e12\)。
这里首先要知道一件事情：之所以一个分数\(\frac{x}{y}\)能够化成有限小数，是因为能够在保证分子为整数的情况下把化成\(10^k\)的形式，也就是说，化成最简分数之后，只有当\(x=y\)或者\(y\)只有素因子\(2,5\)时才能将分数转化为有限小数。
这样的话，朴素的想法是枚举\(x,y\)，然后利用\(\gcd\)化成最简，再进行素因数分解，复杂度\(O(n^2log_2n)\) ，预计得分\(40pts\)。
由于分数化到最简之后分母只可能含有\(2,5\)这两种素因子，所以我们可以考虑先基于优先队列从小到大生成所有的\(y\)，其中\(y<=n\)且\(y\)只可能含有\(2,5\)这两种素因子，显然，如果这些\(y\)作为分母，那么分子是多少都可以化为有限小数。这样的话，我们似乎可以基于最简分数，生成那些不是最简的分数。
在做了一些数学题之后，我知道了一般套路是先写出最后结果的表达式，然后再化简式子简化运算。所以我先给出这个题目的结果\(ans=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[\frac{i}{j}\ is\ finite\ decimal]\) 。感觉先枚举分母，然后枚举合法的分子比较好，所以先交换一下和号：\(ans=\sum\limits_{j=1}^{n}\sum\limits_{i=1}^{n}[\frac{i}{j}\ is\ finite\ decimal]\) 。根据前面的分析，当\(j=1\)或者\(j=2^p5^q\)时，分子随便选，否则，分子就得保证把分母除了\(2,5\)之外的素因子消掉。这给我们一个启示：先枚举分母，然后通过除以\(2,5\)来判断分母的类型，如果是第一种，就直接加结果就好了，如果是第二种，就枚举剩余的真分母的倍数，由于我们的分母\(j\)是不同的，所以虽然可能得到的分数值是一样的，但是分数约化之前仍然是不同的分数。用数学式子表达就是\(ans=\sum\limits_{j=1}^{n}\sum\limits_{i=1}^{[\frac{2^p5^qn}{j}]}1\) ，复杂度接近于$O(n) $ 。预期得分\(80pts\)。
\(subtask\ 3\)的数据范围是\(1e12\)，所以我们可能需要分块。刚才那个式子写得简练一点就是\(ans=\sum\limits_{i=1}^{n}[\frac{2^p5^qn}{i}]\) 。 像这种式子都是可以考虑用分块求和的。但是带着和\(i\)相关的\(2^p5^q\)，不太好处理，我们争取把\(2^p5^q\)拿到前面去。基于这样的愿望，我们可以把式子变成这样：\(ans=\sum2^p5^q\sum\limits_{i=1}^{n}[\frac{n}{i}\ if\ i\%2!=0\ and\ i\%5!=0]\) 。这个式子的意思是，只枚举把分母的\(2,5\)素因子全部提出来的数，然后做乘法。满足这个条件的\(i\)的特点便是\(i\)不含有\(2\)和\(5\)这两种素因子，所以\(i\)的合法的范围并不是\([1,n]\)，看起来又没法直接整除分块。对于每一个合法的\(i\)，我们都要求出来\([\frac{n}{i}]\)，然后再乘上\([1,\lfloor\frac{n}{i}\rfloor]\)之内只含有\(2,5\)这两种因子的数的个数，就是\([\frac{n}{i}]\)的贡献了。只含有\(2,5\)这两种素因子的数可以预处理出来，并且求个数的前缀和。不过由于外层循环还是枚举的\([1,n]\)，所以还是没实质性降低复杂度。
下面想怎么有针对性地枚举\(i\)或者如何分块。我们记\(f(n)\)为\([1,n]\)中的只含有\(2,5\)这两种素因子的数的个数，则\(ans=\sum\limits_{i=1}^{n}f([\frac{n}{i}])[\frac{n}{i}\ if\ i\%2!=0\ and\ i\%5!=0]\) 。我们套路地考虑整除分块，如果全部合法，则\([l,r]\)区间的贡献当然就是\((r-l+1)*f([\frac{n}{i}])*[\frac{n}{i}]\) ，现在由于不合法，所以就不能无脑乘\(r-l+1\)，需要用一些方法统计一下\([l,r]\)中有用的数的个数。什么叫做有用的数呢？就是\([l,r]\)中不含\(2,5\)两种素因子的数的个数。我们知道，区间\([l,r]\)中含有某个素因子\(x\)的数的个数是\([\frac{r}{x}]-[\frac{l-1}{x}]\) （下面那个细节题里面知道的），用所有的数减去这个，就是不含有\(x\)的数的个数了，如果有多个条件，就容斥，总之这不是一件难事儿。注意生成数的时候不要用priority_queue，因为会生成重复的。
代码如下：
#include <bits/stdc++.h> #define ll long long using namespace std; const int N=1e7+9; const ll M=1e12+9; ll only_2_5[N],cnt,n; void get_2_5(); ll getf(ll x); ll getg(ll l,ll r); int main(){ ll ans=0,l=1,r,f,g; scanf("%lld",&n); get_2_5(); for(l=1;l<=n;l=r+1){ r=(n/(n/l)); f=getf(n/l); g=getg(l,r); ans+=f*g*(n/l); } printf("%lld\n",ans); return 0; } ll getg(ll l,ll r){ ll g1,g2,g5,g10; g1=r-l+1; g2=(r/2)-(l-1)/2; g5=(r/5)-(l-1)/5; g10=(r/10)-(l-1)/10; return g1-g2-g5+g10; } ll getf(ll x){ return upper_bound(only_2_5,only_2_5+cnt,x)-only_2_5; //求f(x) } void get_2_5(){ for(ll i=1;i<M;i*=2){ for(ll j=i;j<M;j*=5){ only_2_5[cnt++]=j; } } sort(only_2_5,only_2_5+cnt); }本题收获：
- 整除分块与容斥原理结合时怎么用
- 如何生成\([1,n]\)内所有只含有某些素因子的数
-
题目大意：给定 \(N, M\)，求 \(1 \leq x \leq N\) ,\(1 \leq y \leq M\) 且 \(\gcd(x, y)\) 为质数的 \((x, y)\) 有多少对。
是不是贼熟悉？那就对了，前面有道题和这个描述差不多。显然，这道题的答案是\(ans=\sum\limits_{p\ \in\ prime}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j)=p]\) 。如果暴力计算的话，复杂度是\(O(tn^3logn)\) ，对于\(1e7\)的数据量以及\(1e4\)次询问是完全不可以接受的，属于无脑暴力，我甚至不确定这种做法能拿到分数。这个时候，我们要使用莫比乌斯反演简化计算。同时，这也是我第一次做莫比乌斯反演的题目。
为了使用莫比乌斯反演，我们需要先构造两个数论函数。题解说有一个套路是：设\(f(d)\)为\(\gcd(i,j)=d\)的数对\((i,j)\)的个数，则显然\(f(d)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j)=d]\) 。要计算\(f(d)\)，我们需要二重循环，实在太慢，我们想优化。考虑到我们想用莫比乌斯反演，所以我们必须再构造一个函数，要么新函数求和等于\(f(n)\)，要么\(f(d)\)求和得到新函数。在这里，我们使用了倍数的莫比乌斯反演公式，设\(F(d)\)表示\(\gcd(i,j)=kd,k\in\Z^+\)的数对\((i,j)\)的个数，则\(F(d)=\sum\limits_{d|k\ and\ k\leq\ min(m,n)}{f(k)}\) 。根据乘法原理，可以计算出\(F(d)=[\frac{n}{d}][\frac{m}{d}]\) 。由莫比乌斯反演，知\(f(d)=\sum\limits_{d|k}\mu(\frac{k}{d})F(k)\) ，把\(F(k)\)带入就有\(f(d)=\sum\limits_{d|k}\mu(\frac{k}{d})[\frac{n}{k}][\frac{m}{k}]\) 。所以，$ ans=\sum\limits_{p\ \in\ prime}\sum\limits_{p|k}\mu(\frac{k}{p})[\frac{n}{k}][\frac{m}{k}]$ 。如果计算这个式子，那么复杂度就是\(O(tnlnn)\) ，比刚才好了不少，但还是太慢，能够拿到\(30pts\)，属于高级暴力。
在继续优化之前，我认为我们值得思考为什么通过莫比乌斯反演能够降低复杂度。我们本来是要计算\(f(n)\)的，但是利用莫比乌斯反演，我们利用\(F(n)\)就能表示出\(f(n)\)，而$F(n) $ 的计算甚为简单，只需要\(O(1)\)的时间就好了，所以相对于裸的暴力优化了\(O(n^2)\)的复杂度。
接下来继续考虑优化。对于$ ans=\sum\limits_{p\ \in\ prime}\sum\limits_{p|k}\mu(\frac{k}{p})F(k)$ ，先换一下元，得到\(ans=\sum\limits_{p\ \in\ prime}\sum\limits_{k=1}^{\frac{min(m,n)}{p}}\mu(k)F(kp)\) ，再换元，把\(kp\)换成\(T\)，则\(ans=\sum\limits_{p\ \in\ prime}\sum\limits_{p|T}^{min(n,m)}\mu(\frac{T}{p})F(T)\) ，注意到这个式子里面\(F\)函数计算时只需要考虑\(T\)的取值，所以考虑改变枚举顺序，先枚举\(T\)再枚举\(p\)，这样做的好处是能把\(\mu\)单独放在后面，而\(\mu\)函数是积性函数，有良好的性质，可以简化计算。在改变枚举顺序之后，\(ans=\sum\limits_{T=1}^{min(n,m)}F(T)\sum\limits_{p\ \in\ prime\ and\ p|T}\mu(\frac{T}{p})\) 。 对于内层求和，基于线性筛的思想可以预处理出来，即设\(g(T)=\sum\limits_{p\ \in\ prime\ and\ p|T}\mu(\frac{T}{p})\)，预处理出\(g\)函数的值，这就是\(\mu\)函数的妙处。这样的话，复杂度将降为\(O(tn)\)。但是由于是多组数据访问，所以还是过不了。本做法能够拿到\(60pts\)的好成绩了，属于高级暴力。
如果想继续优化，我们就不能老老实实地枚举。注意到\(F(T)=[\frac{n}{T}][\frac{m}{T}]\) ，可以整除分块，所以，我们对\(g\)函数再求一次前缀和，就可以和\(F\)函数一起利用分块了，这样复杂度\(O(t\sqrt n)\)
代码如下：
#include<bits/stdc++.h> #define ll long long using namespace std; const int N=1e7+9; const int M=1e7+9; ll prime[N],cnt,t,n,m,mobius[N],ans,maxk,f[N],sum[N]; bool isprime[N]; void getprime(); int main(){ scanf("%lld",&t); getprime(); while(t--){ ans=0; scanf("%lld %lld",&n,&m); maxk=min(n,m); // // for(ll i=0;i<cnt;i++){ // for(ll j=1;j*prime[i]<=maxk;j++){ // ans+=mobius[j]*(n/(j*prime[i]))*(m/(j*prime[i])); // } // } // printf("%lld\n",ans); //30pts // for(ll i=1;i<=maxk;i++){ // ans+=(n/i)*(m/i)*f[i]; // } //60pts ll l,r; for(l=1;l<=maxk;l=r+1){ r=min(n/(n/l),(m/(m/l))); ans+=(n/l)*(m/l)*(sum[r]-sum[l-1]); } //100pts printf("%lld\n",ans); } return 0; } void getprime(){ mobius[1]=1; for(ll i=2;i<N;i++){ isprime[i]=true; } for(ll i=2;i<N;i++){ if(isprime[i]){ prime[cnt++]=i; mobius[i]=-1; } for(ll j=0;j<cnt && prime[j]*i<N;j++){ isprime[i*prime[j]]=false; mobius[i*prime[j]]=-mobius[i]; if(i%prime[j]==0){ //某个素因子次数>1 mobius[i*prime[j]]=0; break; } } } for(ll i=0;i<cnt;i++){ for(ll j=1;j*prime[i]<N;j++){ f[j*prime[i]]+=mobius[j]; } } for(ll i=1;i<N;i++){ sum[i]=sum[i-1]+f[i]; } }这个题目给我们的启示有：
- 使用莫比乌斯反演能简化和gcd有关的统计和计算问题，其原理是用一个计算量小的函数去表示计算量大的函数
- 使用莫比乌斯反演时，需要合理地设函数，一种设函数的套路是设\(f(d)\)为\(\gcd(i,j)=d\) 的数对的个数，然后用约数或者倍数的莫比乌斯去构造另一个函数。
- \(\mu\)函数是积性函数，所以可以在线性筛中进行预处理，甚至可以利用线性筛筛除数时的特点做一些奇怪的统计。\(\mu\)函数单身比在一起可能更好用一些。
- 改变求和顺序看似毫无作用，其实可能能够把一些计算式从内层循环拿出到外层循环，这样有利于各自处理各自的和。
- 整除分块是个好东西。
- 如果用莫比乌斯反演不想设函数，可以考虑这个性质\(\sum\limits_{d|n}\mu(d)=[n=1]\)，用\(\gcd(i,j)\)代换\(n\)，有\(\sum\limits_{d|\gcd(i,j)}\mu(d)=[\gcd(i,j)=1]\) 。然后，本题就转换为了\(ans=\sum\limits_{p\ \in \ prime}\sum\limits_{i=1}^{[\frac{n}{p}]}\sum\limits_{j=1}^{[\frac{m}{p}]}[\gcd(i,j)=1]\) ，利用刚才的等式，\(ans=\sum\limits_{p\ \in \ prime}\sum\limits_{i=1}^{[\frac{n}{p}]}\sum\limits_{j=1}^{[\frac{m}{p}]} \sum\limits_{d|\gcd(i,j)}\mu(d)\) 。改变循环顺序，有\(ans=\sum\limits_{p\ \in \ prime}\sum\limits_{d=1}^{min([\frac{n}{p}],[\frac{m}{p}])}\mu(d)\sum\limits_{d|i}\sum\limits_{d|j}1=\sum\limits_{p\ \in \ prime}\sum\limits_{d=1}^{min([\frac{n}{p}],[\frac{m}{p}])}\mu(d)[\frac{n}{d}][\frac{m}{d}]\) 。 最外层\(O(n)\)枚举质数没法优化，但是内部循环可以用分块来优化，只要预处理好莫比乌斯函数的前缀和就好了，总复杂度\(O(n\sqrt n )\)。
- 换元之后再枚举也是一种比较重要的技术，便于分离变量和简化枚举
-
题目大意：对于\(a\leq x\leq b,c\leq y\leq d\)，统计\(\gcd(x,y)=k\)的数对\((x,y)\)的对数。
这个题可以用容斥原理做，并且不涉及什么细节。但是这个题也可以直接去做。我是直接去做的，然后犯了很多错误，下面总结一下。
问题1：在区间\([m,n]\)中，到底有多少个能够被\(k\) 整除的数？
显然，能够被\(k\)整除的最大的数肯定是\(k*[\frac{n}{k}]\)，那最小的呢？是\(k*[\frac{m}{k}]\)吗？显然，\([\frac{m}{k}]\leq\frac{m}{k}\) ，所以\(k*[\frac{m}{k}]\leq m\)，并不一定符合题意。分情况讨论：如果\(k|m\)，则最小的就是它了，那么在这个区间里就有\([\frac{n}{k}]-[\frac{m}{k}]+1\)个这样的数；否则，\(k*[\frac{m}{k}]\)就不是合法的，这时\(k*([\frac{m}{k}]+1)\)是最小的，所以总个数是\([\frac{n}{k}]-[\frac{m}{k}]\) 个。显然我们不想分情况，想把两个结果合成一个。对于\(k|m\)的情况，考虑\([\frac{n}{k}]-[\frac{m-1}{k}]\)，这个时候除非\(k=1\)，否则\(k|(m-1)\)一定不成立，并且，\(k*[\frac{m-1}{k}]\)不在区间内部，是非法的。但是由于计数是头尾都要算，所以\(k|m\)时，结果就是\([\frac{n}{k}]-[\frac{m-1}{k}]\) 。当\(k|m\)不成立时，也考虑那个式子，如果\(k|(m-1)\)也不成立，那么结果也是\([\frac{n}{k}]-[\frac{m-1}{k}]\) ，但若\(k|(m-1)\)成立，那么由于减法只算头尾中的一个，所以也能巧妙地把答案写成\([\frac{n}{k}]-[\frac{m-1}{k}]\) 。综上，有\([\frac{n}{k}]-[\frac{m-1}{k}]\) 个数能被\(k\)整除。
问题2：为什么在整除分块时，块的右端点\(r=min([\frac{n}{[\frac{n}{i}]}],[\frac{m}{[\frac{m}{i}]}])\) ？在这里，我们想求的是类似于\(\sum\limits_{i=1}^{maxi}[\frac{n}{i}][\frac{m}{i}]\) 的式子。
解析：我们之所以想分块，就是想让某段区间内的\([\frac{n}{i}]\)和\([\frac{m}{i}]\)不变，这样才能简化计算。显然，求出的\(r\)端点要同时使两个都在区间内不变，所以要取最小值，这个公式的正确性无可置疑。那么为什么不能先取\(k=min(m,n)\)，然后\(r=[\frac{k}{[\frac{k}{i}]}]\) 呢？考虑\(n=5,m=8,i=3\)，可以计算出\([\frac{n}{[\frac{n}{i}]}]=5,[\frac{m}{[\frac{m}{i}]}]=4\)。所以，必须按照上述的公式计算。
问题3：正数向上取整
推公式时，\(i\)的枚举范围是\([\lceil \frac{a}{k}\rceil,\lfloor\frac{b}{k}\rfloor]\) 。我们当然能手写向上取整的除法，只要以是否整除为分类标准讨论就好了。一定要记住，整数是落在\([a,b]\)区间内的，所以左边向上取整，右边向下取整。
代码如下：
#include <bits/stdc++.h> #define ll long long using namespace std; const int N=5e4+9; const ll o=1; ll prime[N],cnt,mobius[N],sum[N],a,b,c,d,k,n,high,low; bool isprime[N]; void getprime(); int main(){ getprime(); scanf("%lld",&n); while(n--){ scanf("%lld %lld %lld %lld %lld",&a,&b,&c,&d,&k); a=(a-1)/k; b/=k; c=(c-1)/k; d/=k; high=min(b,d); ll ans=0,l=1,r; for(l=1;l<=high;l=r+1){ r=min(b/(b/l),d/(d/l)); if(a >= l) r=min(r,a/(a/l)); //这两行代码我真的不明白怎么回事儿 if(c >= l) r=min(r,c/(c/l)); ans+=(sum[r]-sum[l-1])*(b/l-a/l)*(d/l-c/l); } printf("%lld\n",ans); } return 0; } void getprime(){ mobius[1]=1; for(ll i=2;i<N;i++){ isprime[i]=true; } for(ll i=2;i<N;i++){ if(isprime[i]){ prime[cnt++]=i; mobius[i]=-1; } for(ll j=0;j<cnt && i*prime[j]<N;j++){ isprime[i*prime[j]]=false; mobius[i*prime[j]]=-mobius[i]; if(i%prime[j]==0){ mobius[i*prime[j]]=0; break; } } } for(ll i=1;i<N;i++){ sum[i]=sum[i-1]+mobius[i]; } }另外放上容斥的代码，这个没有任何疑问：
#include <bits/stdc++.h> #define ll long long using namespace std; const int N=5e4+9; const ll o=1; ll prime[N],cnt,mobius[N],sum[N],a,b,c,d,k,n,high,low,ans; bool isprime[N]; void getprime(); ll solve(ll a,ll b); int main(){ getprime(); scanf("%lld",&n); while(n--){ scanf("%lld %lld %lld %lld %lld",&a,&b,&c,&d,&k); ans=0; ans+=solve(b,d); ans-=solve(a-1,d); ans-=solve(b,c-1); ans+=solve(a-1,c-1); printf("%lld\n",ans); } return 0; } ll solve(ll x,ll y){ ll high,ret=0,l=1,r; x=x/k; y=y/k; high=min(x,y); for(l=1;l<=high;l=r+1){ r=min(x/(x/l),y/(y/l)); ret+=(sum[r]-sum[l-1])*(x/l)*(y/l); } return ret; } void getprime(){ mobius[1]=1; for(ll i=2;i<N;i++){ isprime[i]=true; } for(ll i=2;i<N;i++){ if(isprime[i]){ prime[cnt++]=i; mobius[i]=-1; } for(ll j=0;j<cnt && i*prime[j]<N;j++){ isprime[i*prime[j]]=false; mobius[i*prime[j]]=-mobius[i]; if(i%prime[j]==0){ mobius[i*prime[j]]=0; break; } } } for(ll i=1;i<N;i++){ sum[i]=sum[i-1]+mobius[i]; } } -
题目大意：
-
题目大意：给两个超大的整数，求它们的gcd。
先介绍几个算法：
辗转相除法：$\gcd(a,b)=\gcd(b,a%b) $ ，即欧几里得算法，可以用带余除法轻易证明。
更相减损术：\(\gcd(a,b)=\gcd(a,b-a)\)，证明方法和辗转相除法如出一辙，其实辗转相除就是用取模代替了减法。另外，如果\(a,b\)均为偶数，则\(\gcd(a,b)=2\gcd(\frac{a}{2},\frac{b}{2})\) 。
Stein算法：是辗转相除法面对大整数时的改进算法。先说明一个可能用到的结论：若\(\gcd(k,b)=1\)，则\(\gcd(ka,b)=\gcd(a,b)\) ，特殊化之后就是如果一个奇数和一个偶数求\(\gcd\) ，可以先约掉偶数的\(2\)，这样并不影响结果。然后Stein算法的内容是：若\(a,b\)为偶数，则\(\gcd(a,b)=2\gcd(\frac{a}{2},\frac{b}{2})\) ；若\(a,b\)一奇一偶，不妨\(a\)为偶，则\(\gcd(a,b)=\gcd(\frac{a}{2},b)\) ；若\(a,b\)都是奇数，不妨\(a>b\)，则\(\gcd(a,b)=\gcd(\frac{a-b}{2},b)\) ，即更相减损术的内容。显然，每次都至少有一个数大约变成原来的一半，所以整个算法的复杂度也是对数级别的。
这个题我偷懒了，用了python。
# def gcd(a,b): #更相减损术 # while(a!=b): # if(a>b): # a=a-b # else: # b=b-a # return a # def gcd(a,b): #stein，有误 # c=1 # tmp=0 # while(a!=b): # if(a%2==0 and b%2==0): # a=a/2 # b=b/2 # c=c*2 # elif(a%2==0 and b%2!=0): # a=a/2 # elif(a%2!=0 and b%2==0): # b=b/2 # else: # if(a>b): # tmp=b # b=(a-b)/2 # a=tmp # else: # tmp=a # a=(b-a)/2 # b=tmp # return (int)(c*a) def gcd(a,b): #辗转相除法 if(a<b): c=a a=b b=c while(b!=0): a=a%b c=a a=b b=c return a a=(int)(input()) b=(int)(input()) print(gcd(a,b))
组合数学
-
组合数预处理
当规模不太大时，可以考虑\(O(n^2)\)处理
代码如下：
for(int i=0;i<N;i++){ c[i][0]=1; //初始化放在循环里面比较好看 for(int j=1;j<=i;j++){ c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod; //一定注意取模，否则爆ll } }单个组合数计算，可以用定义式，不赘述
处理比较大的组合数，考虑到组合数的定义式，所以我们可以通过处理阶乘的逆元，来做到\(O(1)\) 计算组合数。如P4071 排列计数
显然本题的结果\(ans=C_n^mf(n-m)\) ，其中\(f(n)\)为\(n\)个数的错排方案数，用容斥原理不难计算出\(f(n)=\sum\limits_{i=0}^{n}(-1)^i\frac{n!}{i!}\) 。所以\(ans=\frac{n!}{m!(n-m)!}\sum\limits_{i=0}^{n-m}(-1)^i\frac{(n-m)!}{i!}\)这样的话，我们需要处理\(1!-n!\)的逆元。不过这样的话，单次查询的复杂度还是\(O(n-m)\)，需要继续优化。注意到可以约一下分子分母，得到\(ans=\frac{n!}{m!}\sum\limits_{i=0}^{n-m}(-1)^i\frac{1}{i!}\) ，这样的话，如果预处理出求和部分，就能真正达到\(O(1)\)查询了。
另外出于巩固知识的需要，这里再多说几句：
- 在出现需要对除法的结果取模时（中间取模，防爆ll），我们就必须需要逆元了。
- 线性递推逆元和线性阶乘逆元都要求模数为质数，原因是递推过程中要保证涉及到的数都有逆元，而如果模数是合数，就不一定能保证这一点，即递推时出现了\(inv[0]\)。
- 承接上一条，数论倒数的定义是：如果\(ax\equiv1(mod\ p)\)且\(\gcd(a,p)=1\)，则\(x\)为\(a\)在模\(p\)意义下的逆元。
- 费马小定理求逆元也必须要求模数是质数，这个是由费马小定理本身决定的。
- 扩欧求逆元是不要求\(p\)为质数，但是要求\(\gcd(a,p)=1\)，否则无解。理由的话，参考必要条件：\(ax+by=c\)有整数解的必要条件是\(\gcd(a,b)|c\)，此处\(c=1\)，所以\(\gcd(a,b)\)只能为\(1\) 。
- 由此观之，扩欧求逆元的要求最低，只要存在就可以求；其他三种方法都要求模数为质数才可以求。这个要注意！
关于如何处理阶乘的逆元，首先我们得把\(n!\)算出来，然后用费马小定理或者exgcd干出来，再然后线性倒推回去。线性递推时，考虑到\(inv[i!]=\frac{1}{i!}\)，所以\(inv[(i-1)!]=inv[i!]*i\)
代码如下：
#include <bits/stdc++.h> #define ll long long using namespace std; const int N=1e6+9; const int M=1e6+9; const int mod=1e9+7; ll t,n,m,factinv[N],invsum[N],fact[N]; ll mul(ll a,ll b,ll mod); ll quick_power(ll a,ll b,ll mod); void pre(); int main(){ scanf("%lld",&t); pre(); while(t--){ scanf("%lld%lld",&n,&m); printf("%lld\n",mul(mul(fact[n],factinv[m],mod),invsum[n-m],mod)); } return 0; } void pre(){ ll sign=-1; fact[0]=1; fact[1]=1; for(ll i=2;i<N;i++){ fact[i]=(fact[i-1]*i)%mod; } factinv[N-1]=quick_power(fact[N-1],mod-2,mod); factinv[0]=1; for(ll i=N-1;i>1;i--){ factinv[i-1]=(factinv[i]*i)%mod; } invsum[0]=factinv[0]; for(ll i=1;i<N;i++){ invsum[i]=((invsum[i-1]+sign*factinv[i])%mod+mod)%mod; sign=-sign; } return; } ll quick_power(ll a,ll b,ll mod){ ll base=a,ret=1; while(b){ if(b&1){ ret=mul(ret,base,mod); } base=mul(base,base,mod); b=b>>1; } return ret; } ll mul(ll a,ll b,ll mod){ ll base=a,ret=0; while(b){ if(b&1){ ret=(ret+base)%mod; } base=(base*2)%mod; b=b>>1; } return ret; } -
题目大意：有n个人，m种特产，每种特产有a[i]包，分特产时每个人至少一包，求分配方案数。
我们知道，如果所有特产都是一样的，那么可以转化为不定方程解的组数问题。但是这里有m种不同性质的物品，为了转化，我们先只考虑一种物品。如果先不考虑每个人至少一包，则可以轻易求出来分配某一种特产时的方案数。利用乘法原理，每种特产的方案数乘起来，就是总方案数。下面考虑去除非法情况，枚举有至少1个人，2个人……n个人没有拿到特产的情况，做一个容斥就可以了。
这里我注意到一个细节：由于容斥的时候我们要枚举子集，所以如果没有高效的枚举方案的话，我们要枚举\(2^n\)次（比如带容斥的背包那个题，因为每种硬币不一样贵，所以只能老实枚举），但是这里我们可以依托组合数去枚举，就只需要枚举n次即可了。
-
未完待续
-
这道题太巧妙了，简直是万物皆可容斥！题目大意是：有一个全是0的数表，n行m列，对于每一行每一列，我们都进行且仅进行一次操作：取该行/列的一个前缀，把前缀的所有数加1。最后问我们这个数表一共有多少种情况。
这道题不难发现一个特点：即使是我们给出的操作不同，最终也可能形成相同的数表，其原因在于对列的操作和对行的操作可以“贴边”，比如对第一行的前两个+1配合对第3列第一个+1和对第一行的前三个+1配合第三列前0个加1，得到的就是一样的效果，所以我们可以只算前面那种情况。所以，不是有多少种操作列就有多少种数表。为此我们可以容斥去掉那些重复的。考虑枚举我们的操作序列中至少出现了多少种后面那种情况的操作对，把他们容斥掉就好了。
考虑最后的答案：由于是n行m列，所以每行可以有0,1,2,……m这m+1种操作选择，每一列有n+1种操作选择，所以只考虑操作序列的不同的话，是\((m+1)^n(n+1)^m\)种情况。再把容斥加进去，就是\(ans=\sum\limits_{i=0}^{min(m,n)}(-1)^i(m+1)^{n-i}(n+1)^{m-i}C_n^iC_m^i\)
关于计算，这道题用到了快速幂，且因为数据太大，所以计算组合数用到阶乘逆元。由于涉及取模和加减法，所以最后结果要取模加模再取模。（幸亏模数是质数，所以可以线性做阶乘逆元，出题人真是天使）
代码如下：
#include <bits/stdc++.h> #define ll long long using namespace std; const int N=1e6+9; const int mod=998244353; ll ans=0,fact[N],invfact[N]; void pre(); void solve(); ll quick_power(ll a,ll b,ll mod); int main(){ pre(); solve(); return 0; } void solve(){ ll n,m,k,sign=1,cni,cmi; scanf("%lld %lld",&n,&m); k=min(n,m); for(ll i=0;i<=k;i++){ cni=(((fact[n]*invfact[i])%mod)*invfact[n-i])%mod; cmi=(((fact[m]*invfact[i])%mod)*invfact[m-i])%mod; ans=(ans+sign*(((((fact[i]*quick_power(n+1,m-i,mod)%mod))*quick_power(m+1,n-i,mod)%mod)*cni)%mod*cmi)%mod)%mod; sign=-sign; } ans=((ans%mod)+mod)%mod; printf("%lld\n",ans); } void pre(){ fact[0]=1; fact[1]=1; for(ll i=2;i<N;i++){ fact[i]=(fact[i-1]*i)%mod; } invfact[N-1]=quick_power(fact[N-1],mod-2,mod); for(ll i=N-1;i>0;i--){ invfact[i-1]=(invfact[i]*i)%mod; } } ll quick_power(ll a,ll b,ll mod){ ll base=a,ret=1; while(b){ if(b&1){ ret=(ret*base)%mod; } base=(base*base)%mod; b=b>>1; } return ret; } -
题目大意：给定n，k，找k个不同的组合数\(C_a^b\)，其中\(0\leq b\leq a\leq n\) ，使得这k个组合数不相同，且和最大。
显然，n范围小的时候，可以直接暴力算出所有组合数（二维数组压成一维数组），然后排序，找最后面的k个就可以了。
n范围大的时候应该怎么做呢？通过画杨辉三角形可以发现，最大的组合数是\(C_n^{\frac{n}{2}}\) ，然后它的左右两边是次大的，或许我们可以按照这种策略逐渐把组合数都选出来。这个的思路就是：借助优先队列，首先把a取各种值时最大的组合数塞进去，然后每弹出一个最大的，就把和它a相同但比他小的最大的那个数再塞进去。
这道题有一个难点在于如何比较两个日大的组合数的大小。因为平时算组合数都是取模之后的，所以没法比大小了。在这里我们采用了取对数的方法，把最大的组合数的大小降低到了1e8之内，这样可以用double来存组合数，精度可以接受。在计算组合数时，观察一下表达式发现可以用前缀和优化，这样才能使复杂度正确从而通过此题。
代码如下：
#include <bits/stdc++.h> #define ll long long using namespace std; const int N=1e6+9; const int M=2e6+9; const int mod=1e9+7; typedef struct{ ll a,b,sign; double logvalue; }Node; typedef struct{ ll a,b; }NB; priority_queue<Node> q; set<NB> s; ll n,k,fact[N],invfact[N]; double sum[M],lg[M]; void pre(); void solve(); ll quick_power(ll a,ll b,ll mod); int main(){ scanf("%lld %lld",&n,&k); pre(); solve(); return 0; } void pre(){ fact[0]=1; fact[1]=1; for(int i=2;i<N;i++){ fact[i]=(fact[i-1]*i)%mod; } invfact[N-1]=quick_power(fact[N-1],mod-2,mod); for(ll i=N-1;i>0;i--){ invfact[i-1]=(invfact[i]*i)%mod; } for(ll i=1;i<M;i++){ lg[i]=log(i); } for(ll i=1;i<M;i++){ sum[i]=sum[i-1]+lg[i]; } return ; } ll quick_power(ll a,ll b,ll mod){ ll ret=1,base=a; while(b){ if(b&1){ ret=(ret*base)%mod; } b=b>>1; base=(base*base)%mod; } return ret; } bool operator <(const Node &x,const Node &y){ return x.logvalue<y.logvalue; } void solve(){ Node tmp,now; ll cnt=0,ans=0; for(int i=0;i<=n;i++){ //先把每个a下最大的组合数加入到queue中 now.a=i; now.b=i/2; now.sign=0; now.logvalue=sum[now.a]-sum[now.b]-sum[now.a-now.b]; //前缀和优化 q.push(now); } while(!q.empty()){ if(cnt==k) break; now=q.top(); q.pop(); ans=(ans+((fact[now.a]*invfact[now.b])%mod)*invfact[now.a-now.b])%mod; cnt++; if(now.sign==0){ //描述按照从中间到两边的规则往优先队列里面加数据 if(now.b+1<=now.a){ tmp.a=now.a; tmp.b=now.b+1; tmp.sign=1; tmp.logvalue=sum[tmp.a]-sum[tmp.b]-sum[tmp.a-tmp.b]; q.push(tmp); } } else if(now.sign>0){ tmp.sign=-now.sign; if(now.a/2+tmp.sign>=0){ tmp.a=now.a; tmp.b=now.a/2+tmp.sign; tmp.logvalue=0; tmp.logvalue=sum[tmp.a]-sum[tmp.b]-sum[tmp.a-tmp.b]; q.push(tmp); } } else { tmp.sign=-now.sign+1; if(now.a/2+tmp.sign<=now.a){ tmp.a=now.a; tmp.b=now.a/2+tmp.sign; tmp.logvalue=0; tmp.logvalue=sum[tmp.a]-sum[tmp.b]-sum[tmp.a-tmp.b]; q.push(tmp); } } } printf("%lld\n",ans); } -
题目大意：见原题。
可以发现，我们从右下角开始走，一直走到左上角，所有的这些路径都满足集合的包含关系。
然后我们又注意到一件事情：对于集合\(A\)，如果从中拿出任意一个子集\(B\)，能够快速算出来\(A\)的所有子集中有多少能包含\(B\)就好了。这个是容易计算的：设\(A\) 有n个元素，\(B\)有m个元素，则相当于我们在除了\(B\)中的元素里随便选元素，所以就有\(2^{n-m}\)个符合要求的集合。
先考虑一个比较简单的情况：\(k=2\) 时，就三个集合，我们不难得出这个时候的答案为\(ans=\sum\limits_{i=0}^{n}C_n^i\sum\limits_{j=0}^{n-i}C_{n-i}^{j}\sum\limits_{m=0}^{n-i-j}C_{n-i-j}^{m}\) 。看起来很有意思的一个式子，我不会化简，所以我选择打表观察。（不过后来我推出来了，连续使用二项式定理往回推就行了）打表之后，发现规律是\(2^{2n}\)。于是我们可以猜测，\(n\)个数选\(k\)行的结果是\(2^{nk}\)。结果真的是这样。
-
题目大意：给定n个数对\((a_i,b_i)\)，求\(\sum\limits_{i=0}^{n}\sum\limits_{j=i+1}^{n}C_{a_i+a_j+b_i+b_j}^{a_i+a_j}\) 。\(n\leq 2e5，a_i,b_i\leq 2000\) 。
看到这个问题，直接计算肯定不行。要想简化运算，必须对式子进行变形。考虑里面那个组合数的实际意义：有一个方程为\(x_1+x_2+……+x_{a_i+a_j+1}=b_i+b_j\) ，求其非负解的组数。额好像也没看出来什么。要想加速，可能需要考虑一下子计算出好几个组合数的值。
然后不会。
看了题解，发现确实是考虑实际意义，只不过是往几何的方面想的。考虑\(C_{a+b}^{a}\) ,其意义是，从\((0,0)\)走到\((a,b)\)，有多少种走法。在这里就是从\((0,0)\)到\((a_i+a_j,b_i+b_j)\) 有多少种走法。再转化一步，就是\((a_i,b_i)\)和\((-a_j,-b_j)\) 之间的路径的条数。直接计算还是单独算组合数，不能加快，得考虑几何性质。
可以发现，从第三象限的点到第一象限的点，其路径必然会穿过\(y=-x\)这条直线，并且每条路径仅通过其中一个点。所以一个计算方法是：统计经过\(y=-x\)上的每个点的路径有多少条，求和即可。这样的复杂度就降到了\(O(n)\) 。
如何统计呢？考虑经过某个叫做\((a,-a)\)的点的路径的条数，显然，如果是从\((-a_j,-b_j)\)走到\((a_i,b_i)\)的话，只要\(-a_j\leq a\leq a_i ,-b_j\leq-a\leq b_i\) ，就是可以走到的。有多少条呢？两个组合数一算就有了。
-
组合数奇偶性结论：对于\(C_n^k\)，如果\(n\&k=k\)，则组合数为奇数。
线性代数
-
三根柱子的汉诺塔，只能相邻柱移动，有n个盘子，求移动步数。
n个盘子要想移动完，需要先把n-1个盘子从第一根柱子移到第三个柱子，然后把第n个盘子移到第二个柱子，再把n-1个盘子移到第一个柱子，再把第n个盘子移到第三个柱子，再把n-1个盘子移到第三个柱子，所以，\(f[n]=3f[n-1]+2\) 。这个题n的规模是\(1e9\)，所以要加速转移过程。两个方式：一是求出通项公式，然后快速幂，二是利用矩阵快速幂优化dp。
-
主要是代码，先放一个低精度版本的代码：
#include <bits/stdc++.h> #define ll long long #define eps 0.00001 using namespace std; const int N=1e2+9; double a[N][N]; double b[N]; double ans[N]; ll n; bool flag; int main(){ scanf("%lld",&n); for(int i=1;i<=n;i++){ for(int j=1;j<=n;j++){ scanf("%lf",&a[i][j]); } scanf("%lf",&b[i]); } for(int i=1;i<=n;i++){ if(fabs(a[i][i])<eps){ //如果系数是0 flag=false; for(int j=i+1;j<=n;j++){ if(fabs(a[j][i])>eps){ swap(a[i],a[j]); swap(b[i],b[j]); flag=true; break; } } } else { //如果系数不是0 flag=true; } if(flag){ for(int j=i+1;j<=n;j++){ //对下面各行进行消元 if(fabs(a[j][i])>=eps){ //如果j行需要消元 double coef=a[j][i]/a[i][i]; //求出要乘的系数 for(int k=1;k<=n;k++){ a[j][k]=a[j][k]-a[i][k]*coef; //消元 } b[j]=b[j]-b[i]*coef; //增广部分消元 } } } } //至此消元完成 for(int i=n;i>=1;i--){ double ret=b[i]; for(int j=i+1;j<=n;j++){ ret=ret-a[i][j]*ans[j]; } if(fabs(a[i][i])<eps){ printf("No Solution\n"); return 0; } ans[i]=ret/a[i][i]; } for(int i=1;i<=n;i++){ printf("%.2lf\n",ans[i]); } return 0; } -
算法大家都懂，一个原矩阵一个单位阵做相同操作，最后原来的单位阵被化成的矩阵就是逆矩阵了。看看代码就好：
#include <bits/stdc++.h> #define ll long long using namespace std; const int mod=1e9+7; const int N=404; typedef struct{ ll matrix[N][N]; ll n; void elmination(ll x,ll y,ll k){ for(ll i=1;i<=n;i++){ matrix[y][i]=((matrix[y][i]-matrix[x][i]*k)%mod+mod)%mod; } } void permutation(ll x,ll y){ swap(matrix[x],matrix[y]); } void mul(ll x,ll k){ for(ll i=1;i<=n;i++){ matrix[x][i]=(matrix[x][i]*k)%mod; } } }Matrix; Matrix a,b; ll n; ll quick_power(ll a,ll b); int main(){ scanf("%lld",&n); a.n=n; b.n=a.n; for(int i=1;i<=n;i++){ for(int j=1;j<=n;j++){ scanf("%lld",&a.matrix[i][j]); } b.matrix[i][i]=1; } for(int i=1;i<=n;i++){ //化为上三角形 bool flag=true; if(a.matrix[i][i]==0){ flag=false; for(int j=i+1;j<=n;j++){ if(a.matrix[j][i]!=0){ a.permutation(i,j); b.permutation(i,j); flag=true; break; } } } else { flag=true; } if(flag==false){ printf("No Solution\n"); return 0; } ll inv=quick_power(a.matrix[i][i],mod-2); for(int j=i+1;j<=n;j++){ ll tmp=a.matrix[j][i]; a.elmination(i,j,(inv*tmp)%mod); b.elmination(i,j,(inv*tmp)%mod); } } // for(int i=1;i<=n;i++){ // for(int j=1;j<=n;j++){ // printf("%lld ",b.matrix[i][j]); // } // printf("\n"); // } for(int i=n;i>=1;i--){ ll inv=quick_power(a.matrix[i][i],mod-2); for(int j=i-1;j>=1;j--){ ll tmp=a.matrix[j][i]; //注意这里很容易错！！！ a.elmination(i,j,(inv*tmp)%mod); b.elmination(i,j,(inv*tmp)%mod); } a.mul(i,inv); b.mul(i,inv); } for(int i=1;i<=n;i++){ for(int j=1;j<=n;j++){ printf("%lld ",b.matrix[i][j]); } printf("\n"); } return 0; } ll quick_power(ll a,ll b){ ll ret=1,base=a; while(b){ if(b&1){ ret=(ret*base)%mod; } base=(base*base)%mod; b=b>>1; } return ret; } -
又是一个模板题。这个考察了线性方程组解的情况的判断。这道题还比较良心，给的是\(n\)个\(n\)元一次方程，不需要再讨论秩与未知数个数的问题了。
假设系数矩阵是\(n\)行\(m\)列的，其秩为\(r\)，方程等号右边组成向量\(b\)，对于一般的情形，我们分别讨论一下解的情况
-
若\(r=m<n\) ，说明未知数个数小于方程个数，可能会出现矛盾方程。这时候，系数矩阵化简成对角形之后的形状如下所示（我们假设\(r=2\)进行演示）：
\[\left( \begin{matrix} 1 & 0 \\ 0 & 1 \\ … & … \\ … & … \end{matrix} \right) \]缺省部分全部都用\(0\)来填补，因为只要有前两行在，后面的肯定都能被第一行的多少倍与第二行的多少倍消成\(0\) 。
因为下面全都是\(0\)了，所以如果要有解，就要对\(b\)有要求了。如果\(b\)经过同样变换之后，下面出现了非\(0\)分量，从方程上看就是出现了矛盾方程，那么方程组无解，否则，方程组有解，且仅有唯一解。
-
若\(r=n<m\)，说明方程个数小于未知数个数，每一行都有一个主元，这样化成对角形矩阵的时候形状类似于下面所示
\[\left( \begin{matrix} 1 & 0 & …\\ 0 & 1 & … \end{matrix} \right) \]缺省部分在这里就不一定全是\(0\)了。
这种情况下，显然后面的那些自由变量是可以随便取值的，不管如何取值我们都可以解出来非自由变量的值，所以这种情况下有无穷多组解。
这个时候我们发现自由变量的系数们都是\(0\)了，也就是说自由变量任意地取值对于非自由变量的解来说都没有影响，而显然非自由变量肯定能解出一组解来。所以，这种情况下方程组有无穷多组解。
-
若\(r<min(n,m)\) ，则矩阵化简出来类似这样：
\[\left( \begin{matrix} 1 & 0 & …\\ 0 & 1 & …\\ … & … & … \end{matrix} \right) \]缺省部分不一定全是\(0\)，但是保证前两列下面的缺省部分全是\(0\)。
显然这个时候不一定有解，根据前面的分析，是否有解要看下面的\(n-r\)行中的\(b\)向量在变换之后有没有非\(0\)元，有的话，那就无解了；没有的话，保证有解，但解的个数还不能确定。在保证有解的情况下，根据前面的分析，自由变量可以随便取值，所以有无穷多组解。
-
若\(r=n=m\) ，矩阵化为最简形就是单位阵。根据前面的分析，对于任意的\(b\)，都有解；由于没有自由变量，所以解唯一。
综上所述：
-
\(r=m<n\) 时：
- \(b\)向量选择得当，则有唯一解。
- \(b\)向量选择不得当，则无解。
-
\(r=n<m\)时：
方程组有无穷多解。
-
\(r<min(n,m)\)时：
- \(b\)向量选择得当，则有无穷多解。
- \(b\)向量选择不得当，则无解。
-
\(r=n=m\)时：
必定有唯一解。
容易发现，我们在叙述中提到了\(b\)向量选择“得当”这样一个模糊的概念，经过观察，我们发现，\(b\)向量选择得当，意味着解方程组时行变换之后没有矛盾方程，在矩阵形状上看就是系数矩阵右边再加上一列\(b\)之后，与只有系数矩阵相比，在进行同样的初等行变换之后不出现系数矩阵某行全为0，而b向量变换后在该行不为0的情况，即系数矩阵的秩\(=\)增广矩阵的秩。因此，\(b\)向量选择得当\(<=>\)系数矩阵的秩\(=\)增广矩阵的秩。
所以我们还能准确地得出一个关于线性方程组是否有解的结论：
- 系数矩阵的秩\(=\)增广矩阵的秩时，方程组一定有解。
- 否则，方程组一定无解。
