数学
高精度（大数）¶
- 注意L（数组长度会极大的影响算法效率，应根据需要的数据范围设置）
+¶
const int L=11010;
string add(string a,string b)//只限两个非负整数相加
{
string ans;
int na[L]={0};
int nb[L]={0};
int la=a.size();
int lb=b.size();
for(int i=0;i<la;i++) na[la-1-i]=a[i]-'0';
for(int i=0;i<lb;i++) nb[lb-1-i]=b[i]-'0';
int lmax=max(la,lb);
for(int i=0;i<lmax;i++) na[i]+=nb[i],na[i+1]+=na[i]/10,na[i]%=10;
if(na[lmax]) lmax++;
for(int i=lmax-1;i>=0;i--) ans+=na[i]+'0';
return ans;
}
-¶
const int L=10010;
string sub(string a,string b)//只限大的非负整数减小的非负整数
{
string ans;
int na[L]={0};
int nb[L]={0};
int la=a.size();
int lb=b.size();
for(int i=0;i<la;i++) na[la-1-i]=a[i]-'0';
for(int i=0;i<lb;i++) nb[lb-1-i]=b[i]-'0';
int lmax=max(la,lb);
for(int i=0;i<lmax;i++)
{
na[i]-=nb[i];
if(na[i]<0) na[i]+=10,na[i+1]--;
}
while(!na[--lmax]&&lmax>0) lmax++;
for(int i=lmax-1;i>=0;i--) ans+=na[i]+'0';
return ans;
}
*¶
const int L=100005;
string mul(string a,int b)//高精度a乘单精度b
{
int na[L];
string ans;
int la=a.size();
fill(na,na+L,0);
for(int i=la-1;i>=0;i--) na[la-i-1]=a[i]-'0';
int w=0;
for(int i=0;i<la;i++) na[i]=na[i]*b+w,w=na[i]/10,na[i]=na[i]%10;
while(w) na[la++]=w%10,w/=10;
la--;
while(la>=0) ans+=na[la--]+'0';
return ans;
}
const int L=110;
string mul(string a,string b)
{
string s;
int na[L],nb[L],nc[L],La=a.size(),Lb=b.size();
fill(na,na+L,0);fill(nb,nb+L,0);fill(nc,nc+L,0);
for(int i=La-1;i>=0;i--) na[La-i]=a[i]-'0';
for(int i=Lb-1;i>=0;i--) nb[Lb-i]=b[i]-'0';
for(int i=1;i<=La;i++)
for(int j=1;j<=Lb;j++)
nc[i+j-1]+=na[i]*nb[j];
for(int i=1;i<=La+Lb;i++)
nc[i+1]+=nc[i]/10,nc[i]%=10;
if(nc[La+Lb]) s+=nc[La+Lb]+'0';
for(int i=La+Lb-1;i>=1;i--)
s+=nc[i]+'0';
return s;
}
- FTT优化版（\(O(n^2)\)->\(O(nlog(n))\)）
#define L(x) (1 << (x))
const double PI = acos(-1.0);
const int Maxn = 133015;
double ax[Maxn], ay[Maxn], bx[Maxn], by[Maxn];
char sa[Maxn/2],sb[Maxn/2];
int sum[Maxn];
int x1[Maxn],x2[Maxn];
int revv(int x, int bits)
{
int ret = 0;
for (int i = 0; i < bits; i++)
{
ret <<= 1;
ret |= x & 1;
x >>= 1;
}
return ret;
}
void fft(double * a, double * b, int n, bool rev)
{
int bits = 0;
while (1 << bits < n) ++bits;
for (int i = 0; i < n; i++)
{
int j = revv(i, bits);
if (i < j)
swap(a[i], a[j]), swap(b[i], b[j]);
}
for (int len = 2; len <= n; len <<= 1)
{
int half = len >> 1;
double wmx = cos(2 * PI / len), wmy = sin(2 * PI / len);
if (rev) wmy = -wmy;
for (int i = 0; i < n; i += len)
{
double wx = 1, wy = 0;
for (int j = 0; j < half; j++)
{
double cx = a[i + j], cy = b[i + j];
double dx = a[i + j + half], dy = b[i + j + half];
double ex = dx * wx - dy * wy, ey = dx * wy + dy * wx;
a[i + j] = cx + ex, b[i + j] = cy + ey;
a[i + j + half] = cx - ex, b[i + j + half] = cy - ey;
double wnx = wx * wmx - wy * wmy, wny = wx * wmy + wy * wmx;
wx = wnx, wy = wny;
}
}
}
if (rev)
{
for (int i = 0; i < n; i++)
a[i] /= n, b[i] /= n;
}
}
int solve(int a[],int na,int b[],int nb,int ans[])
{
int len = max(na, nb), ln;
for(ln=0; L(ln)<len; ++ln);
len=L(++ln);
for (int i = 0; i < len ; ++i)
{
if (i >= na) ax[i] = 0, ay[i] =0;
else ax[i] = a[i], ay[i] = 0;
}
fft(ax, ay, len, 0);
for (int i = 0; i < len; ++i)
{
if (i >= nb) bx[i] = 0, by[i] = 0;
else bx[i] = b[i], by[i] = 0;
}
fft(bx, by, len, 0);
for (int i = 0; i < len; ++i)
{
double cx = ax[i] * bx[i] - ay[i] * by[i];
double cy = ax[i] * by[i] + ay[i] * bx[i];
ax[i] = cx, ay[i] = cy;
}
fft(ax, ay, len, 1);
for (int i = 0; i < len; ++i)
ans[i] = (int)(ax[i] + 0.5);
return len;
}
string mul(string sa,string sb)
{
int l1,l2,l;
int i;
string ans;
memset(sum, 0, sizeof(sum));
l1 = sa.size();
l2 = sb.size();
for(i = 0; i < l1; i++)
x1[i] = sa[l1 - i - 1]-'0';
for(i = 0; i < l2; i++)
x2[i] = sb[l2-i-1]-'0';
l = solve(x1, l1, x2, l2, sum);
for(i = 0; i<l || sum[i] >= 10; i++) // 进位
{
sum[i + 1] += sum[i] / 10;
sum[i] %= 10;
}
l = i;
while(sum[l] <= 0 && l>0) l--; // 检索最高位
for(i = l; i >= 0; i--) ans+=sum[i] + '0'; // 倒序输出
return ans;
}
/¶
using namespace std;
string div(string a,int b)//高精度a除以单精度b
{
string r,ans;
int d=0;
if(a=="0") return a;//特判
for(int i=0;i<a.size();i++)
{
r+=(d*10+a[i]-'0')/b+'0'; //求出商
d=(d*10+(a[i]-'0'))%b; //求出余数
}
int p=0;
for(int i=0;i<r.size();i++)
if(r[i]!='0') {p=i;break;}
return r.substr(p);
}
const int L=110;
int sub(int *a,int *b,int La,int Lb)
{
if(La<Lb) return -1;//如果a小于b，则返回-1
if(La==Lb)
{
for(int i=La-1;i>=0;i--)
if(a[i]>b[i]) break;
else if(a[i]<b[i]) return -1;//如果a小于b，则返回-1
}
for(int i=0;i<La;i++)//高精度减法
{
a[i]-=b[i];
if(a[i]<0) a[i]+=10,a[i+1]--;
}
for(int i=La-1;i>=0;i--)
if(a[i]) return i+1;//返回差的位数
return 0;//返回差的位数
}
string div(string n1,string n2,int nn)//n1,n2是字符串表示的被除数，除数,nn是选择返回商还是余数
{
string s,v;//s存商,v存余数
int a[L],b[L],r[L],La=n1.size(),Lb=n2.size(),i,tp=La;//a，b是整形数组表示被除数，除数，tp保存被除数的长度
fill(a,a+L,0);fill(b,b+L,0);fill(r,r+L,0);//数组元素都置为0
for(i=La-1;i>=0;i--) a[La-1-i]=n1[i]-'0';
for(i=Lb-1;i>=0;i--) b[Lb-1-i]=n2[i]-'0';
if(La<Lb || (La==Lb && n1<n2)) {
//cout<<0<<endl;
return n1;}//如果a<b,则商为0，余数为被除数
int t=La-Lb;//除被数和除数的位数之差
for(int i=La-1;i>=0;i--)//将除数扩大10^t倍
if(i>=t) b[i]=b[i-t];
else b[i]=0;
Lb=La;
for(int j=0;j<=t;j++)
{
int temp;
while((temp=sub(a,b+j,La,Lb-j))>=0)//如果被除数比除数大继续减
{
La=temp;
r[t-j]++;
}
}
for(i=0;i<L-10;i++) r[i+1]+=r[i]/10,r[i]%=10;//统一处理进位
while(!r[i]) i--;//将整形数组表示的商转化成字符串表示的
while(i>=0) s+=r[i--]+'0';
//cout<<s<<endl;
i=tp;
while(!a[i]) i--;//将整形数组表示的余数转化成字符串表示的</span>
while(i>=0) v+=a[i--]+'0';
if(v.empty()) v="0";
//cout<<v<<endl;
if(nn==1) return s;
if(nn==2) return v;
}
max¶
string max(const string& a, const string& b) {
if (a.length() > b.length()) {
return a;
} else if (a.length() < b.length()) {
return b;
}
for (std::size_t i = 0; i < a.length(); ++i) {
if (a[i] > b[i]) {
return a;
} else if (a[i] < b[i]) {
return b;
}
}
return a;
}
min¶
string min(const std::string& a, const std::string& b) {
if (a.length() < b.length()) {
return a;
} else if (a.length() > b.length()) {
return b;
}
for (std::size_t i = 0; i < a.length(); ++i) {
if (a[i] < b[i]) {
return a;
} else if (a[i] > b[i]) {
return b;
}
}
return a;
}
mod¶
int mod(string a,int MOD)
{
int left=0;
for(int i=0;i<a.size();i++)
left=(left*10+(a[i]-'0'))%MOD;
return left;
}
int main()
{
string a; int b;
while(cin>>a>>b) cout<<mod(a,b)<<endl;
return 0;
}
阶乘¶
const int L=100005;
int a[L];
string fac(int n)
{
string ans;
if(n==0) return "1";
fill(a,a+L,0);
int s=0,m=n;
while(m) a[++s]=m%10,m/=10;
for(int i=n-1;i>=2;i--)
{
int w=0;
for(int j=1;j<=s;j++) a[j]=a[j]*i+w,w=a[j]/10,a[j]=a[j]%10;
while(w) a[++s]=w%10,w/=10;
}
while(!a[s]) s--;
while(s>=1) ans+=a[s--]+'0';
return ans;
}
幂¶
#define L(x) (1 << (x))
const double PI = acos(-1.0);
const int Maxn = 133015;
double ax[Maxn], ay[Maxn], bx[Maxn], by[Maxn];
char sa[Maxn/2],sb[Maxn/2];
int sum[Maxn];
int x1[Maxn],x2[Maxn];
int revv(int x, int bits)
{
int ret = 0;
for (int i = 0; i < bits; i++)
{
ret <<= 1;
ret |= x & 1;
x >>= 1;
}
return ret;
}
void fft(double * a, double * b, int n, bool rev)
{
int bits = 0;
while (1 << bits < n) ++bits;
for (int i = 0; i < n; i++)
{
int j = revv(i, bits);
if (i < j)
swap(a[i], a[j]), swap(b[i], b[j]);
}
for (int len = 2; len <= n; len <<= 1)
{
int half = len >> 1;
double wmx = cos(2 * PI / len), wmy = sin(2 * PI / len);
if (rev) wmy = -wmy;
for (int i = 0; i < n; i += len)
{
double wx = 1, wy = 0;
for (int j = 0; j < half; j++)
{
double cx = a[i + j], cy = b[i + j];
double dx = a[i + j + half], dy = b[i + j + half];
double ex = dx * wx - dy * wy, ey = dx * wy + dy * wx;
a[i + j] = cx + ex, b[i + j] = cy + ey;
a[i + j + half] = cx - ex, b[i + j + half] = cy - ey;
double wnx = wx * wmx - wy * wmy, wny = wx * wmy + wy * wmx;
wx = wnx, wy = wny;
}
}
}
if (rev)
{
for (int i = 0; i < n; i++)
a[i] /= n, b[i] /= n;
}
}
int solve(int a[],int na,int b[],int nb,int ans[])
{
int len = max(na, nb), ln;
for(ln=0; L(ln)<len; ++ln);
len=L(++ln);
for (int i = 0; i < len ; ++i)
{
if (i >= na) ax[i] = 0, ay[i] =0;
else ax[i] = a[i], ay[i] = 0;
}
fft(ax, ay, len, 0);
for (int i = 0; i < len; ++i)
{
if (i >= nb) bx[i] = 0, by[i] = 0;
else bx[i] = b[i], by[i] = 0;
}
fft(bx, by, len, 0);
for (int i = 0; i < len; ++i)
{
double cx = ax[i] * bx[i] - ay[i] * by[i];
double cy = ax[i] * by[i] + ay[i] * bx[i];
ax[i] = cx, ay[i] = cy;
}
fft(ax, ay, len, 1);
for (int i = 0; i < len; ++i)
ans[i] = (int)(ax[i] + 0.5);
return len;
}
string mul(string sa,string sb)
{
int l1,l2,l;
int i;
string ans;
memset(sum, 0, sizeof(sum));
l1 = sa.size();
l2 = sb.size();
for(i = 0; i < l1; i++)
x1[i] = sa[l1 - i - 1]-'0';
for(i = 0; i < l2; i++)
x2[i] = sb[l2-i-1]-'0';
l = solve(x1, l1, x2, l2, sum);
for(i = 0; i<l || sum[i] >= 10; i++) // 进位
{
sum[i + 1] += sum[i] / 10;
sum[i] %= 10;
}
l = i;
while(sum[l] <= 0 && l>0) l--; // 检索最高位
for(i = l; i >= 0; i--) ans+=sum[i] + '0'; // 倒序输出
return ans;
}
string Pow(string a,int n)
{
if(n==1) return a;
if(n&1) return mul(Pow(a,n-1),a);
string ans=Pow(a,n/2);
return mul(ans,ans);
}
gcd¶
const int L=110;
string add(string a,string b)
{
string ans;
int na[L]={0},nb[L]={0};
int la=a.size(),lb=b.size();
for(int i=0;i<la;i++) na[la-1-i]=a[i]-'0';
for(int i=0;i<lb;i++) nb[lb-1-i]=b[i]-'0';
int lmax=la>lb?la:lb;
for(int i=0;i<lmax;i++) na[i]+=nb[i],na[i+1]+=na[i]/10,na[i]%=10;
if(na[lmax]) lmax++;
for(int i=lmax-1;i>=0;i--) ans+=na[i]+'0';
return ans;
}
string mul(string a,string b)
{
string s;
int na[L],nb[L],nc[L],La=a.size(),Lb=b.size();//na存储被乘数，nb存储乘数，nc存储积
fill(na,na+L,0);fill(nb,nb+L,0);fill(nc,nc+L,0);//将na,nb,nc都置为0
for(int i=La-1;i>=0;i--) na[La-i]=a[i]-'0';//将字符串表示的大整形数转成i整形数组表示的大整形数
for(int i=Lb-1;i>=0;i--) nb[Lb-i]=b[i]-'0';
for(int i=1;i<=La;i++)
for(int j=1;j<=Lb;j++)
nc[i+j-1]+=na[i]*nb[j];//a的第i位乘以b的第j位为积的第i+j-1位（先不考虑进位）
for(int i=1;i<=La+Lb;i++)
nc[i+1]+=nc[i]/10,nc[i]%=10;//统一处理进位
if(nc[La+Lb]) s+=nc[La+Lb]+'0';//判断第i+j位上的数字是不是0
for(int i=La+Lb-1;i>=1;i--)
s+=nc[i]+'0';//将整形数组转成字符串
return s;
}
int sub(int *a,int *b,int La,int Lb)
{
if(La<Lb) return -1;//如果a小于b，则返回-1
if(La==Lb)
{
for(int i=La-1;i>=0;i--)
if(a[i]>b[i]) break;
else if(a[i]<b[i]) return -1;//如果a小于b，则返回-1
}
for(int i=0;i<La;i++)//高精度减法
{
a[i]-=b[i];
if(a[i]<0) a[i]+=10,a[i+1]--;
}
for(int i=La-1;i>=0;i--)
if(a[i]) return i+1;//返回差的位数
return 0;//返回差的位数
}
string div(string n1,string n2,int nn)//n1,n2是字符串表示的被除数，除数,nn是选择返回商还是余数
{
string s,v;//s存商,v存余数
int a[L],b[L],r[L],La=n1.size(),Lb=n2.size(),i,tp=La;//a，b是整形数组表示被除数，除数，tp保存被除数的长度
fill(a,a+L,0);fill(b,b+L,0);fill(r,r+L,0);//数组元素都置为0
for(i=La-1;i>=0;i--) a[La-1-i]=n1[i]-'0';
for(i=Lb-1;i>=0;i--) b[Lb-1-i]=n2[i]-'0';
if(La<Lb || (La==Lb && n1<n2)) {
//cout<<0<<endl;
return n1;}//如果a<b,则商为0，余数为被除数
int t=La-Lb;//除被数和除数的位数之差
for(int i=La-1;i>=0;i--)//将除数扩大10^t倍
if(i>=t) b[i]=b[i-t];
else b[i]=0;
Lb=La;
for(int j=0;j<=t;j++)
{
int temp;
while((temp=sub(a,b+j,La,Lb-j))>=0)//如果被除数比除数大继续减
{
La=temp;
r[t-j]++;
}
}
for(i=0;i<L-10;i++) r[i+1]+=r[i]/10,r[i]%=10;//统一处理进位
while(!r[i]) i--;//将整形数组表示的商转化成字符串表示的
while(i>=0) s+=r[i--]+'0';
//cout<<s<<endl;
i=tp;
while(!a[i]) i--;//将整形数组表示的余数转化成字符串表示的</span>
while(i>=0) v+=a[i--]+'0';
if(v.empty()) v="0";
//cout<<v<<endl;
if(nn==1) return s;
if(nn==2) return v;
}
bool judge(string s)//判断s是否为全0串
{
for(int i=0;i<s.size();i++)
if(s[i]!='0') return false;
return true;
}
string gcd(string a,string b)//求最大公约数
{
string t;
while(!judge(b))//如果余数不为0，继续除
{
t=a;//保存被除数的值
a=b;//用除数替换被除数
b=div(t,b,2);//用余数替换除数
}
return a;
}
压位完整高精¶
#include <algorithm>
using namespace std;
const int N_huge = 850, base = 100000000;
struct huge {
typedef long long value;
value a[N_huge];
int len;
void clear() { len = 1; a[len] = 0; }
huge() { clear(); }
huge(value x) { *this = x; }
huge operator =(huge b) {
len = b.len;
for (int i = 1; i <= len; ++i)
a[i] = b.a[i];
return *this;
}
huge operator +(huge b) {
int L = len > b.len ? len : b.len;
huge tmp;
for (int i = 1; i <= L + 1; ++i)
tmp.a[i] = 0;
for (int i = 1; i <= L; ++i) {
if (i > len)
tmp.a[i] += b.a[i];
else if (i > b.len)
tmp.a[i] += a[i];
else {
tmp.a[i] += a[i] + b.a[i];
if (tmp.a[i] >= base) {
tmp.a[i] -= base;
++tmp.a[i + 1];
}
}
}
if (tmp.a[L + 1])
tmp.len = L + 1;
else
tmp.len = L;
return tmp;
}
huge operator -(huge b) {
int L = len > b.len ? len : b.len;
huge tmp;
for (int i = 1; i <= L + 1; ++i)
tmp.a[i] = 0;
for (int i = 1; i <= L; ++i) {
if (i > b.len)
b.a[i] = 0;
tmp.a[i] += a[i] - b.a[i];
if (tmp.a[i] < 0) {
tmp.a[i] += base;
--tmp.a[i + 1];
}
}
while (L > 1 && !tmp.a[L])
--L;
tmp.len = L;
return tmp;
}
huge operator *(huge b) {
int L = len + b.len;
huge tmp;
for (int i = 1; i <= L; ++i)
tmp.a[i] = 0;
for (int i = 1; i <= len; ++i)
for (int j = 1; j <= b.len; ++j) {
tmp.a[i + j - 1] += a[i] * b.a[j];
if (tmp.a[i + j - 1] >= base) {
tmp.a[i + j] += tmp.a[i + j - 1] / base;
tmp.a[i + j - 1] %= base;
}
}
tmp.len = len + b.len;
while (tmp.len > 1 && !tmp.a[tmp.len])
--tmp.len;
return tmp;
}
pair<huge, huge> divide(huge a, huge b) {
int L = a.len;
huge c, d;
for (int i = L; i; --i) {
c.a[i] = 0;
d = d * base;
d.a[1] = a.a[i];
int l = 0, r = base - 1, mid;
while (l < r) {
mid = (l + r + 1) >> 1;
if (b * mid <= d)
l = mid;
else
r = mid - 1;
}
c.a[i] = l;
d -= b * l;
}
while (L > 1 && !c.a[L])
--L;
c.len = L;
return make_pair(c, d);
}
huge operator /(value x) {
value d = 0;
huge tmp;
for (int i = len; i; --i) {
d = d * base + a[i];
tmp.a[i] = d / x;
d %= x;
}
tmp.len = len;
while (tmp.len > 1 && !tmp.a[tmp.len])
--tmp.len;
return tmp;
}
value operator %(value x) {
value d = 0;
for (int i = len; i; --i)
d = (d * base + a[i]) % x;
return d;
}
huge operator /(huge b) { return divide(*this, b).first; }
huge operator %(huge b) { return divide(*this, b).second; }
huge &operator +=(huge b) { *this = *this + b; return *this; }
huge &operator -=(huge b) { *this = *this - b; return *this; }
huge &operator *=(huge b) { *this = *this * b; return *this; }
huge &operator ++() {
huge T;
T = 1;
*this = *this + T;
return *this;
}
huge &operator --() {
huge T;
T = 1;
*this = *this - T;
return *this;
}
huge operator ++(int) {
huge T, tmp = *this;
T = 1;
*this = *this + T;
return tmp;
}
huge operator --(int) {
huge T, tmp = *this;
T = 1;
*this = *this - T;
return tmp;
}
huge operator +(value x) {
huge T;
T = x;
return *this + T;
}
huge operator -(value x) {
huge T;
T = x;
return *this - T;
}
huge operator *(value x) {
huge T;
T = x;
return *this * T;
}
huge operator *=(value x) { *this = *this * x; return *this; }
huge operator +=(value x) { *this = *this + x; return *this; }
huge operator -=(value x) { *this = *this - x; return *this; }
huge operator /=(value x) { *this = *this / x; return *this; }
huge operator %=(value x) { *this = *this % x; return *this; }
bool operator ==(value x) {
huge T;
T = x;
return *this == T;
}
bool operator !=(value x) {
huge T;
T = x;
return *this != T;
}
bool operator <=(value x) {
huge T;
T = x;
return *this <= T;
}
bool operator >=(value x) {
huge T;
T = x;
return *this >= T;
}
bool operator <(value x) {
huge T;
T = x;
return *this < T;
}
bool operator >(value x) {
huge T;
T = x;
return *this > T;
}
huge operator =(value x) {
len = 0;
while (x) {
a[++len] = x % base;
x /= base;
}
if (!len)
a[++len] = 0;
return *this;
}
bool operator <(huge b) {
if (len < b.len)
return 1;
if (len > b.len)
return 0;
for (int i = len; i; --i) {
if (a[i] < b.a[i])
return 1;
if (a[i] > b.a[i])
return 0;
}
return 0;
}
bool operator ==(huge b) {
if (len != b.len)
return 0;
for (int i = len; i; --i)
if (a[i] != b.a[i])
return 0;
return 1;
}
bool operator !=(huge b) { return !(*this == b); }
bool operator >(huge b) { return !(*this < b || *this == b); }
bool operator <=(huge b) { return (*this < b) || (*this == b); }
bool operator >=(huge b) { return (*this > b) || (*this == b); }
void str(char s[]) {
int l = strlen(s);
value x = 0, y = 1;
len = 0;
for (int i = l - 1; i >= 0; --i) {
x = x + (s[i] - '0') * y;
y *= 10;
if (y == base)
a[++len] = x, x = 0, y = 1;
}
if (!len || x)
a[++len] = x;
}
void read() { scanf("%s", s); this->str(s); }
void print() {
printf("%d", (int)a[len]);
for (int i = len - 1; i; --i) {
for (int j = base / 10; j >= 10; j /= 10) {
if (a[i] < j)
printf("0");
else
break;
}
printf("%d", (int)a[i]);
}
printf("\n");
}
};
例题¶
组合数学¶
- 排列
- \(P_n^r=\frac{n!}{(n-r)!}\)
- 圆排列\(\frac{P^r_n}{r}\)
- 组合
- \(C_n^r=\frac{n!}{r!(n-r)!}\)
- 性质公式
- \(C_n^r=C_n^{n-r}\)
- \(C_n^r=C_{n-1}^r+C_{n-1}^{r-1}\)帕斯卡公式，为了避免阶层及逆元时使用dp思想快速计算组合数
- \(C_n^0+\dots+C_n^n=2^n\)
- 多重集排列
- k个不同元素，每个有无数个，选r个排列\(k^r\)
- k个不同元素，分别\(n_1,\dots,n_k\)个（共n个），排列个数为\(\frac{n!}{n_1!\dots n_k!}\)
二项式定理（杨辉三角）¶
- 用于递推计算组合数（帕斯卡公式）
- 二项式定理\((a+b)^n=\sum_{r=0}^nC_n^ra^rb^{n-r}\)
计算组合数¶
-
使用递推计算
-
计算多项，复杂度\(O(n^2)\)
-
```c++ const int N = 1005; #define mod 10007 int c[N][N]; int dfs(int n,int m){ //用递推公式求二项式系数 if(!m) return c[n][m]=true; if(m==1) return c[n][m]=n; if(c[n][m]) return c[n][m]; //记忆化 if(n-m<m) m=n-m; return c[n][m]=(dfs(n-1,m)+dfs(n-1,m-1))%mod;
}for (int i = 0; i < N; i++) c[i][0] = 1; for (int i = 1; i < N; i++) for (int j = 1; j <= i; j++) c[i][j] = (c[i - 1][j] + c[i - 1][j - 1])%mod;
```
-
用逆直接计算
-
只计算1项，复杂度\(O(n)\)
-
```c++ #define mod 10007 int fac[10001]; //预计算阶乘 int inv[10001]; //预计算逆 int fastPow(int a, int n){} //和上一个代码完全一样 int C(int n,int m){ //算组合数，用到除法取模的逆 return (fac[n]inv[m]%modinv[n-m]%mod)%mod; } int main(){ int a,b,n,m,k,ans; cin >>n>>r; fac[0] = 1; for(int i=1;i<=n;i++){ fac[i] = (fac[i-1]*i)%mod; //预计算阶乘，要取模 inv[i] = fastPow(fac[i],mod-2); //用费马小定理预计算i!的逆 } ans = C(n,r); cout << ans; return 0; }
```
卢卡斯定理¶
-
由于逆元可能无法求出，递推求解组合数又太慢，卢卡斯定理用于特定情况快速求解组合数
-
如果m小于n可能在求解\(r \ n-r\)得逆元时无法满足互素得条件
-
或者卢卡斯定理也可以用于解决求解超大组合数的问题（\(O(n)\)也超时）
-
要求m是一个素数
-
\(C_n^r\mod m=\prod_{i=0}^kC_{n_i}^{r_i}\mod m\)
-
\(n=n_km^k+\dots+n_1m+n_0 \ r=r_km^k+\dots+r_1m+r_0\)是m进制展开
-
如果存在\(n_i<r_i\)有\(c_{n_i}^{r_i}=0\)即\(c_n^i \mod m=0\)
-
更常用得表述\(C_n^r\equiv C_{n\mod m}^{r\mod m}C_{n/m}^{r/m}(\mod m)\)递归求解
-
由于\(r\mod m\)和\(n \mod m\)显然小于m，这两部分可以直接求逆计算
-
复杂度约为\(log_m^n\)
const int N = 100010;
typedef long long ll;
ll fac[N]; //预计算阶乘，取模
ll fastPow(ll a, ll n, ll m){ //标准快速幂
ll ans = 1;
a %= m; //防止下面的ans*a越界
while(n) {
if(n & 1) ans = (ans*a) % m;
a = (a*a) % m;
n >>= 1;
}
return ans;
}
ll inverse(ll a,int m){ return fastPow(fac[a],m-2,m); } //用费马小定理计算逆
ll C(ll n,ll r,int m){ //用逆计算C(n mod m, r mod m) mod m
if(r>n)return 0;
return ((fac[n] * inverse(r,m))%m * inverse(n-r,m)%m);
}
ll Lucas(ll n,ll r,int m){ //用递归计算C(n, r) mod m
if(r==0) return 1;
return C(n%m,r%m,m) * Lucas(n/m,r/m,m)%m; //分两部分计算
}
int main(){
int T; cin >> T;
while(T--){
int a,b,m; cin >>a>>b>>m;
fac[0] = 1;
for(int i=1;i<=m;i++) fac[i]=(fac[i-1]*i)%m; //预计算阶乘，取模
cout << Lucas(a+b,a,m) << endl;
}
}
```
- 如果m不是一个素数，可以尝试将m进行质因数分解，然后使用中国剩余定理求解
- <img src="https://thdlrt.oss-cn-beijing.aliyuncs.com/image-20230923104243801.png" alt="image-20230923104243801" style="zoom:33%;" />等价求解
- <img src="https://thdlrt.oss-cn-beijing.aliyuncs.com/image-20230923104305186.png" alt="image-20230923104305186" style="zoom:33%;" />
### 鸽巢原理
- 把$kn+1$个物体放到n个盒子，至少有1个盒子至少包含$k+1$个物体
- $q_1+q_2+\dots+q_n-n+1$个物体放入n个盒子，那么或者第一个至少有$q_1$个物体，要么第二个盒子至少有$1_2\dots$
- Ramsey定理：6人（或更多）中，或者三个人两两认识，或者三个人都不认识。
#### 例子
- n个人，认识的人两两握手，至少有两个人握手次数相同
- K种糖果，每个给定数目，判断能不能不两天吃一样的
- 数目最多得糖果作为隔板
- 5个自然数，必有三个数和被3整除
- 相似问题可以计算前缀和，按照前缀和模得余数分类，找相同类别得两点作为区间
### 容斥原理
- $|A_1\bigcup A_2\bigcup\dots\bigcup A_n|=\sum|A_i|-\sum|A_i\bigcap A_j|+\dots+(-1)^{n+1}|A_1\bigcap\dots\bigcap A_n|$
#### 优化dp
- 可以使用二进制表示状态（状态压缩），枚举差分
- 用全部方案数目减去不合法方案数目得到何方方案数目
- [P1450 硬币购物 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)](https://www.luogu.com.cn/problem/P1450)
- 容斥原理优化背包dp
- 将困难得多重背包转化为一系列完全背包得组合
- 而可以通过预处理完全背包，快速完成每次查询
- [P5664 Emiya 家今天的饭 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)](https://www.luogu.com.cn/problem/solution/P5664)
### 卡特兰数
- $C_n=\frac{1}{n+1}C_{2n}^n \ n=0,1,2.\dots$
- $1,1,2,5,14,42,132,429,1430\dots$
#### 求解
- $C_n=C_{2n}^n-C_{2n}^{n-1}$
- $O(n)$
- 预处理阶乘、逆元，快速计算
- ```c++
ll fac[2000001], inv[2000001],mod=20100403;
ll fastPow(ll n,ll num)
{
ll ans=1,sum=n;
while(num)
{
if(num&1)
{
ans=(ans*sum)%mod;
}
sum=(sum*sum)%mod;
num>>=1;
}
return ans;
}
int C(int n,int m){
return (fac[n]*inv[m]%mod*inv[n-m]%mod)%mod;
}
int main() {
ll n;
cin>>n;
fac[0] = 1;
for(int i=1;i<=(n+m);i++){
fac[i] = (fac[i-1]*i)%mod;
inv[i] = fastPow(fac[i],mod-2);
}
ll res=((C(n*2,n)-C(n*2,n-1))%mod+mod)%mod;
cout<<res<<endl;
return 0;
}
```
- 递推$C_n=\sum C_kC_{n-k} \ \ \ C_0=1$
- $O(n^2)$
- 分治模型
- 无需除法，编码简单
- ```c++
ll C[n+1];
memset(C,0,sizeof(C));
C[0]=C[1]=1;
for(int i=2;i<=n;i++){
for(int j=0;j<i;j++){
C[i]=(C[i]+(C[j]*C[i-j])%mod)%mod;
}
}
```
- 递推$C_n=\frac{4n-2}{n+1}C_{n-1} \ C_0=1$
- $O(n)$
- 可以看出卡特兰数增长非常快$O(4^n)$
- 方法1、3相对块很多，但需要预处理阶层、计算逆元等
#### 应用
- n\*n 棋盘在不穿过主对角线（始终在右下方）的情况下的走法
- 跨过对角线的路径可以映射为一种新的路径
- <img src="https://thdlrt.oss-cn-beijing.aliyuncs.com/image-20230922214306874.png" alt="image-20230922214306874" style="zoom:33%;" />
- 2n个位置放n个1、0，要求任何前缀都是1得数目大于等于0（前缀和，把0看作-1，1看作1）
- <img src="https://thdlrt.oss-cn-beijing.aliyuncs.com/image-20230922213233205.png" alt="image-20230922213233205" style="zoom: 33%;" />
- 可以扩展为括号匹配问题、数入栈出栈产生的序列问题（出栈操作小于等于入栈）
- 也可以将0、1转化为向右和向上走，转化为棋盘问题
- <img src="https://thdlrt.oss-cn-beijing.aliyuncs.com/image-20230922214908220.png" alt="image-20230922214908220" style="zoom: 25%;" />
- 一个矩形只能覆盖阶梯的一个角，因此左下角的矩形也一定覆盖一个，由此可以将问题拆分为两个子问题，得到类似$C_n=\sum C_kC_{n-k}$
- n个节点构成的不同二叉树的数目
- 以根节点进行划分，也可以得到类似$C_n=\sum C_kC_{n-k}$
- n+1条边的凸多边形区域内插入不相交的边，划分为三角形的方案数$C_{n-1}$
### 斯特林数
- 为指数增长，通常需要配合高精度运算
- 区别：圆排列中有顺序而盒子没有
#### 第一类斯特林数
- $s(n,k)$把n个不同元素分到k个**圆排列**中（不能为空），有多少方案
- $s(n,k)=s(n-1,k-1)+(n-1)s(n-1,k), \ 1<=k<=n$
- $s(0,0)=1,s(k,0)=0, \ 1<=k<=n$
#### 第二类斯特林数
- $S(n,k)$把n个不同的球分配到k个无差别的盒子，不能有空，分配方式数目
- $S(n,k)=S(n-1,k-1)+kS(n-1,k), \ 1<=k<=n$
- $S(0,0)=1,S(i,0)=0, \ 1<=i<=n$
```c++
ll S[N][N];
ll getS(ll n,ll m){
if(S[n][m])
return S[n][m];
if(n==0&&m==0) return 1;
if(m==0||n==0) return 0;
S[n][m]=getS(n-1,m-1)+m*getS(n-1,m);
return S[n][m];
}
- 高精度版本
const int N=105;
const int L=100005;
string mul(string a,int b)
{
int na[L];
string ans;
int la=a.size();
fill(na,na+L,0);
for(int i=la-1;i>=0;i--) na[la-i-1]=a[i]-'0';
int w=0;
for(int i=0;i<la;i++) na[i]=na[i]*b+w,w=na[i]/10,na[i]=na[i]%10;
while(w) na[la++]=w%10,w/=10;
la--;
while(la>=0) ans+=na[la--]+'0';
return ans;
}
string add(string a,string b)//只限两个非负整数相加
{
string ans;
int na[L]={0};
int nb[L]={0};
int la=a.size();
int lb=b.size();
for(int i=0;i<la;i++) na[la-1-i]=a[i]-'0';
for(int i=0;i<lb;i++) nb[lb-1-i]=b[i]-'0';
int lmax=max(la,lb);
for(int i=0;i<lmax;i++) na[i]+=nb[i],na[i+1]+=na[i]/10,na[i]%=10;
if(na[lmax]) lmax++;
for(int i=lmax-1;i>=0;i--) ans+=na[i]+'0';
return ans;
}
string S[N][N];
string getS(int n,int m){
if(S[n][m]!="")
return S[n][m];
if(n==0&&m==0) return "1";
if(m==0||n==0) return "0";
S[n][m]=add(getS(n-1,m-1),mul(getS(n-1,m),m));
return S[n][m];
}
例题¶
- Problem - E - Codeforces
- 对于长度为 \(t\)（不难发现最长路径约为\(\log{n}\)） 的任何路径，城市最大值不大于 \(k\) 的赋值数量是 \(k^t\)（即枚举路径上的最大值）。因此，城市最大值恰好为 \(k\) 的赋值数量是 \(k^t - (k-1)^t\)，而所有赋值中城市最大值的总和是 \(\sum_{k=1}^m(k^t - (k-1)^t)k\)（即\(m^t\times m-\sum_{k=1}^{m-1}k^t\)）。我们可以在 \(\mathcal{O}(m\log n)\) 的时间内解决这个问题。
- 接下来统计长度为t的路径数目（不能使用dfs，点数目为\(O(10^{18})\)）
- 不难发现，本题成城市本质上是一个二叉树（并且除了最后一行外是一个满二叉树），因此每行只有一个节点为根的子树不是完全二叉树，也就是说子树形状的种类的数目是\(O\log{n}\)
- 让 \(dp_{i,j}\) 表示具有 \(i\) 个顶点的子树中长度为 \(j\) 的路径数量，让 \(f_{i,j}\) 表示在具有 \(i\) 个顶点的子树中以 \(i\) 为一端且长度为 \(j\) 的路径数量。(可以统计一个子树从跟出发的各路径长度，由此合并得到个长度路径数目（路径长度也就\(O\log(n)\)，每层枚举组合，\(O\log(n)^3\)）)
数论¶
模运算¶
- 注意模运算可能改变两个数的大小关系
-
有时答案要取正数（尤其是大数减小数）可以使用
(x%mod+mod)%mod取正数 -
负数的模：c++按照正整数求余，符号与被除数保持一致
-
对于两个大数（ll）相乘，及时使用模运算也不能直接相乘否则可能会造成溢出，转化为加法并使用类似快速幂的方式加速(ab->a*2b/2+a*2*b%2)
-
c++ ll mul(ll a,ll b,ll m){ a = a%m; b = b%m; ll res =0; while(b>0){ if(b&1)res=(res+a)%m; a=(a+a)%m; b>>=1; } return res; }
快速幂¶
long long quickpow(int n,long long num)
{
long long ans=1,sum=n;
while(num)
{
if(num&1)
{
ans=(ans*sum)%mod;
}
sum=(sum*sum)%mod;
num>>=1;
}
return ans;
}
龟速乘¶
-
防止两个较大的longlong在乘法过程中溢出
-
c++ ll mul(ll a,ll b,ll m){ //乘法取模：a*b % m(防溢出龟速乘法) ll res=0; while(b>0){ if(b&1) res=(res+a) % m; a=(a+a) % m; b>>=1; } return res; }
gcd与lcm¶
- a能整除b（b%a==0）记作a|b
gcd¶
- 性质
- \(gcd(a,b)=gcd(a,k*a+b)\)
- \(k*gcd(a,b)=gcd(k*a,k*b)\)
- \(gcd(a,b)=d\)则\(gcd(a/d,b/d)=1\)（即互素）
- \(gcd(a,b,c)=gcd(a,gcd(b,c))\)
- \(gcd(a,b)=gcd(a/c,b)\)（\(gcd(b,c)=1\)）
- gcd区间可以使用st表等结构维护
- 如 \(gcd(l,r)=gcd(gcd(l,k_2),gcd(k_1,r))\) 只需要直接对重叠区间求 gcd
- \(gcd(p^{c_{1}},p^{c_{2}})=p^{gcd(c_{1},c_{2})}\)
算法¶
- 欧几里得算法
-
\(O(\log{n})\)n为较大值
ll gcd(ll a,ll b) { return b>0 ? gcd(b,a%b):a; } -
更相损益术 ^95660a
- 避免了大整数取模的性能问题，但是其依靠两数求差的方式来递归，运算次数大于辗转相除法
- 最坏时间复杂度 \(O(n)\)
-
stein算法
- 只需要使用加减法和移位，高精时代码较少
int gcd(int a,int b){ if(a<b) { //arrange so that a>b int temp = a; a = b; b=temp; } if(0==b)//the base case return a; if(a%2==0 && b%2 ==0)//a and b are even return 2*gcd(a/2,b/2); if ( a%2 == 0)// only a is even return gcd(a/2,b); if ( b%2==0 )// only b is even return gcd(a,b/2); return gcd((a+b)/2,(a-b)/2);// a and b are odd } -
\(O(log(max(a, b)))\)
-
最大比率 -求解 (等比数列）分数的最大公约数
- 设比率为 \(\left( \frac{u}{d} \right)^k\)
- 那么待求公因数的项可以表示为 \(\left( \frac{u}{d} \right)^{c_{i}}\)，\(c_{i}\) 是 \(k\) 的倍数，目标就是求出 \(gcd(c_{1}\dots c_{i})\) ，结果就是 \(\left( \frac{u}{d} \right)^{gcd}\)
- 使用更相损益术
- \(p^{gcd(c_{1,c_{2}})}=p^{gcd(c_{1},c_{2}-c_{1})}=gcd\left( p^{c_{1}},\frac{p^{c_{2}}}{p^{c_{1}}} \right)\)
- 直到相除后结果为 1 ：\(gcd(p^{c_{1}},1)=gcd(p^{c_{1}},p^0)=p^{gcd(c_{1},0)}=p^{c_{1}}\)
lcm¶
-
算数基本定理：\(n=p_1^{c_1}p_2^{c_2}...\)pi为素数
-
c++ ll lcm(ll a,ll b){ return a/gcd(a,b)*b; } -
裴蜀定理¶
-
若 ab 均为整数，则有整数 xy 使得 \(ax+by=gcd(a,b)\) 即 \(ax+by\) 是 \(gcd(a,b)\) 的整数倍
-
对于 \(a_{1}x_{1}+\dots+a_{n}x_{n}=c\) 只存在有限个 \(c\) 不存在等式的条件是 \(gcd(a_{1},\dots,a_{n})=1\)
例题¶
-
给出xyk求xy第k大公约数
-
所有公因数都能整除最大公因数
-
因此先找到最大公因数，然后遍历\([1,\sqrt{gcd(x,y)}]\)判断整除即可得到全部公因数
-
hdu4497
-
\(gcd(x,y,z)=G\) \(lcm(x,y,z)=L\) 求满足条件的xyz有多少种
-
由
显然一定要有L%G=0 -
问题转化为\(gcd(x,y,z)=1\) \(lcm(x,y,z)=L/G\)
-
对于每一项质数的分解的幂，要求三个数种至少一个为0（gcd=1），并且最大值必须为L/G对应的值
-
\[ \{0,0,t_1\}3种\\ \{0,t_1,t_1\}3种\\ \{0,t_1,[1,t_1-1]6(t-1)种\} \]
-
即对于该位一共有6*t1种（L/G的质因数分解）
-
-
质因数分解，按照大小关系分类讨论
-
求gcd值一定小于等于原先的a和b，即多所有数求一次gcd即可，注意把负数变正
-
poj1597
- 给定 a b x 要求\((ax)\%b\)可以取得\([0,b-1]\)的全部值
- \(ax+by=d\)，只需要\(d=1\)，即判断\(gcd(a,b)=1\)
线性丢番图方程¶
- \(ax+by=c\) abc为已知整数，求是否存在整数解xy，即直线上是否存在整数坐标点
- 使用裴蜀定理可以知道：如果\(gcd(a,b)=d|c\)则方程有整数解
- 且存在无数多个整数解，可以表示为\(x=x_0+(b/d)n\) \(y=y_0-(a/d)n\)
- 如果要求解数据范围内的解，可以先求出一个整数解，然后派生出范围内的其他解
- 比如求解\(p_1(x_1,y_1)p_2(x_2,y_2)\)线段上两点之间整数点的数目
- 不难求出两点直线为\((y_2-y_1)x+(x_1-x_2)y=y_2x_1-y_1x_2\)
- \(a=y_2-y_1 \\ b=x_1-x_2 \\ c=y_2x_1-y_1x_2\)
- 由此转化为线性丢番图问题
线性丢番图与扩展欧几里得算法¶
-
需要先知道一个格点才能派生出其他点，而求点本质上就是求解\(ax+by=gcd(a,b)\)，可以使用扩展欧几里得算法求解
-
c++ ll extend_gcd(ll a,ll b,ll &x,ll &y){ if(b == 0){ x=1; y=0; return a;} ll d = extend_gcd(b,a%b,y,x); y -= a/b * x; return d; } -
计算gcd的同时记录使得等式成立的\(x_0 \ y_0\)
-
\(ax_0c/d+by_0c/d=d*c/d\)
-
即\(x_o^
=x_0c/d \ y_0^=y_0c/d\) -
\((n-m)t+kL=x-y\)建模
-
c++ #include<bits/stdc++.h> using namespace std; #define ll long long ll extend_gcd(ll a,ll b,ll &x,ll &y){ if(b == 0){ x=1; y=0; return a;} ll d = extend_gcd(b,a%b,y,x); y -= a/b * x; return d; } int main(){ ll n,m,x,y,L; cin>>x>>y>>m>>n>>L; ll a=n-m,c=x-y; if(a<0){ a=-a; c=-c;} //处理负数相当于将两只青蛙的先后对调 ll d = extend_gcd(a,L,x,y); if(c%d != 0) cout<<"Impossible"; //判断方程有无解。 else cout<<((x*(c/d))%(L/d)+(L/d))%(L/d); //x的最小整数解 }- 关于负数：扩展欧几里得求解\(ax+by=c\)时应当保证ab为非负数（否则结果错误），不过要注意的是，如果在计算前取反了，要在得到结果后对相应的x\y取反
同余问题¶
-
若\(m|(a-b)\)称ab对模m同余
-
即\(a=b+km\)
-
完全剩余系：\(modm\)下的完全剩余系为\(1,2,\dots,m-1\)
-
性质：对于\(a\equiv b(modm) \ c\equiv d(modm)\)
-
\[ a+c\equiv b+d(modm)\\ a-c\equiv b-d(modm)\\ ac\equiv bd(modm)\\ a^k \equiv b^k(modm)(k>0) \]
乘法逆元¶
- \(ax\equiv1(modm)\)的解\(x\)为a模m的逆\(a^{-1}\)
- 要求\(gcd(a,m)=1\)（否则无解）
扩展欧几里得¶
-
乘法逆元的求解本质上就是解线性丢番图方程，自然可以用扩展欧几里得解决
-
\(O(logn)\)
-
c++ ll extend_gcd(ll a,ll b,ll &x,ll &y){ if(b == 0){ x=1; y=0; return a;} ll d = extend_gcd(b,a%b,y,x); y -= a/b * x; return d; } long long mod_inverse(long long a, long long m){ long long x,y; extend_gcd(a,m,x,y); return (x%m + m) % m; } //返回1表示无解
欧拉定理¶
-
常用于对指数求模进行化简
-
\(gcd(a,m)=1, \ a^{\varphi(m)}\equiv1(\mod m)\)
-
推论：
-
\[ a^b \equiv \begin{cases} a^{b\mod\varphi(m)} & \text{gcd}(a,m)=1 \\ a^b & b<\varphi(m) \\ a^{b\mod{\varphi(m)}+\varphi(m)} & b\geq\varphi(m) \end{cases} \ (\mod m) \]
费马小定理¶
-
\(a^{m-1}=1(modm)\)m是素数且a与m互素
-
是欧拉定理的一个特例
- 则\(a^{m-2}\)就是\(a\)的逆
-
特殊情况下一种更方便的求法
-
\(O(logm)\)
c++ long lnog mod_inverse(long long a, long long mod){ return fast_pow(a,mod-2,mod); }
递推求解¶
-
求1~n内全部元素的逆
-
设\(p/i=k\)即\(ki+r\equiv0(modp)\)
-
左右乘以\(i^{-1}r^{-1}\)得到\(kr^{-1}+i^{-1}\equiv0(modp)\)
-
这里的-1是指模p下的逆，由于前面的元素已经计算过了，因此可以进行递推
-
\(i^{-1}\equiv-kr^{-1}(modp)\)即\(i^{-1}\equiv (p-p/i)r^{-1}(modp)\)
-
\(O(n)\)
-
c++ long long inv[N]; void inverse(long long n,long long p){ inv[1]=1; for(int i=2;i<N;i++) inv[i]=(p-p/i)*inv[p%i]%p; }
除法取模¶
- \((a/b)modm=(ab^{-1}modm)=(a\mod m)(b^{-1}\mod m)\)
同余方程¶
- 求\(x\)使得\(ax\equiv b(mod m)\)
- 其实就是\(ax-b=km\)等价于线性丢番图方程
- \(gcd(a,m)=d|b\)时有\(d\)个在\(m\)完全剩余系下不同的解\(x=x_0+(m/d)*n(0<=n<d)\)
-
当\(gcd(a,m)=1\)互素时有唯一的解
-
使用逆元求解
- \(x\equiv a^{-1}b(modm)\)
同余方程组¶
-
\[ x\equiv a_1\mod m_1\\ x\equiv a_2\mod m_2\\ \vdots\\ x\equiv a_r\mod m_r \]
-
结果为无穷多\(x=x_0+nm_1\dots m_r\)
中国剩余定理¶
-
要求\(m_1,m_2,\dots,m_r\)两两互素
-
\(M_i=M/m_i \ \ \ M_i^{-1}为M_i模m_i的逆元\)
-
\(M=m_1m_2\dots m_r\)
-
\(x\equiv(a_1M_1M_1^{-1}+\dots +a_rM_rM_r^{-1})(\mod M)\)
-
c++ ll n,a[16],m[16],Mi[16],mul=1,X; void exgcd(ll a,ll b,ll &x,ll &y){ if(b==0){x=1;y=0;return ;} exgcd(b,a%b,x,y); ll z=x;x=y,y=z-y*(a/b); } int main(){ n=rd(); for(int t=1;t<=n;++t){ int M=rd();m[t]=M; mul*=M; a[t]=rd(); } for(int t=1;t<=n;++t){ Mi[t]=mul/m[t]; ll x=0,y=0; exgcd(Mi[t],m[t],x,y); //a[t]=(a[t]%m[t]+m[t])%m[t]; //如果存在a[t]<0 X+=a[t]*Mi[t]*(x<0?x+m[t]:x);//此除可以取%mul防止溢出，也可能需要使用龟速乘 } printf("%lld",X%mul); return 0; }
迭代法¶
-
用于解决m之间不互质的问题，每次合并两个同余式，逐渐合并
-
-
```c++ #include
using namespace std; typedef long long ll; const int N = 100010; int n; ll ai[N], mi[N]; ll mul(ll a,ll b,ll m){ //乘法取模：ab % m(防溢出龟速乘法) ll res=0; while(b>0){ if(b&1) res=(res+a) % m; a=(a+a) % m; b>>=1; } return res; } ll extend_gcd(ll a,ll b,ll &x,ll &y){ //扩展欧几里得 if(b == 0){ x=1; y=0; return a;} ll d = extend_gcd(b,a%b,y,x); y -= a/b * x; return d; } ll excrt(){ //求解同余方程组，返回最小正整数解 ll x,y; ll m1 = mi[1], a1 = ai[1]; //第1个等式 ll ans = 0; for(int i=2;i<=n;i++){ //合并每2个等式
ll a2 = ai[i], m2 = mi[i]; // 第2个等式 //合并为：aX + bY = c
ll a = m1, b = m2, c = (a2 - a1%m2 + m2) % m2; //下面求解 aX + bY = c ll d = extend_gcd(a,b,x,y); //用扩展欧几里得求x0 if(c%d != 0) return -1; //无解 x = mul(x,c/d,b/d); //aX + bY = c 的特解t，最小值
ans = a1 + x m1; //代回原第1个等式，求得特解x' m1 = m2/d*m1; //先除再乘，避免越界。合并后的新m1（即lcm） ans = (ans%m1 + m1) % m1; //最小正整数解 a1 = ans; //合并后的新a1 } return ans; } int main(){ scanf("%d", &n); for(int i=1;i<=n;++i) scanf("%lld%lld",&mi[i],&ai[i]); printf("%lld",excrt()); return 0; }
```
例题¶
-
hdu1576：已知\(n=A\%9973 \ A|B \ \gcd(B,9973)=1\)求解\((A/B)\%9973\)
-
这道题不能使用除法逆元直接求解，需要变形后使用扩展欧几里得解答
-
\[ 设k=(a/b)\%mod\\ A/B=k+9973x\\ 即A=kB+9973xB\\ 带入n得 \ n=KB\%9973\\ kB=n+8873y\\ (k/n)B+(-y/n)9973=1 \ 转化为线性丢番图方程\\ 求解得到k=(k/n)*k \]
-
```c++ int main(void) { int t, i; long n, b, a=9973, x, y;
scanf("%d", &t); for(i=0; i<t; i++) { scanf("%ld%ld", &n, &b); exgcd(b, a, &x, &y); x = (x + a) % a; // x有可能为负，需要调整 printf("%ld\n", x*n%a); } return 0;} ```
- \(n=ax+by(x,y∈N)\)的数是合法的，不能表示的数就是最小能表示的数-1
- 由于\(gcd(a,b)=1\)有解，由此我们也就可以得到\(ax+by=1\)的通解，我们找到这样两个特殊的解\(ax
+by=1\)和\(ax+by=1\)用于与\(n=ax+by(x,y∈N)\)做差，其中\(x`\)是满足的最小非负整数\(y``\)同样。 - 显然\(y` \ x``\)一定是负数，这两项不可能导致小于零
- 我们要关心的是\(a(x-x`)\)和\(b(y-y``)\)只有这两项都小于0才说明无解
- 我们要找的就是\(n=a(x`-1)+b(y``-1) \ ans=n-1\)
- 此外这两项实际上是相邻的，即\(x
=x`-b \ y=y`+a\)，原式可以化简\(ab-a-b\)
素数问题¶
- 伯特兰-切比雪夫定理
- 若整数n > 3，则至少存在一个质数p，符合n < p < 2n − 2。
素数判定¶
-
试除法（适用于较小的数字）
-
复杂度\(O(\sqrt{n})\)
-
c++ bool is_prime(ll n){ if(n<=1) return false; for(ll i=2;i<=sqrt(n);i++) if(n%i==0) return false; return true; } -
Miller-Rabin素性测试
-
基于费马素性测试
- 尝试逆用费马小定理，随机选a如果n满足式子，那么n大概率就是素数，不过在Carmichael（561 1105...）时会出错，此时式子恒成立
-
复杂度\(O(s(log_2n)^3)\)（s通常取值为50）
-
c++ #include <bits/stdc++.h> typedef long long LL; LL fast_pow(LL x,LL y,int m){ //快速幂取模：x^y mod m LL res = 1; x %= m; while(y) { if(y&1) res = (res*x) % m; x = (x*x) % m; y>>=1; } return res; } bool witness(LL a, LL n){ // Miller-Rabin素性测试。返回true表示n是合数 LL u = n-1; //注意，u的意义是：n-1的二进制去掉末尾0 int t = 0; // n-1的二进制，是奇数u的二进制，后面加t个零 while(u&1 == 0) u = u>>1, t++; // 整数n-1末尾0的个数，就是t LL x1, x2; x1 = fast_pow(a,u,n); // 先计算 a^u mod n for(int i=1; i<=t; i++) { // 做t次平方取模 x2 = fast_pow(x1,2,n); // x1^2 mod n if(x2 == 1 && x1 != 1 && x1 != n-1) return true; //用推论判断 x1 = x2; } if(x1 != 1) return true; //用费马测试判断是否为合数：an-1≡1(mod n)不成立，是合数 return false; } int miller_rabin(LL n,int s){ //对n做s次测试 if(n<2) return 0; if(n==2) return 1; //2是素数 if(n % 2 == 0 ) return 0; //偶数 for(int i = 0;i < s && i < n;i++){ //做s次测试 LL a = rand() % (n - 1) + 1; //基值a是随机数 if(witness(a,n)) return 0; //n是合数，返回0 } return 1; //n是素数，返回1 } int main(){ int m; while(scanf("%d",&m) != EOF){ int cnt = 0; for(int i = 0; i < m; i++){ LL n; scanf("%lld",&n); int s = 50; //做s次测试 cnt += miller_rabin(n,s); } printf("%d\n",cnt); } return 0; }
素数筛选¶
-
埃氏筛
-
\(O(nloglogn)\)
```c++ const int N = 1e7; //定义空间大小，1e7约10M int prime[N+1]; //存放素数，它记录visit[i] = false的项 bool visit[N+1]; int E_sieve(int n) { for(int i = 0; i <= n; i++) visit[i]= false; for(int i = 2; ii <= n; i++) //筛掉非素数。改为i<=sqrt(n)，计算更快 if(!visit[i]) for(int j=ii; j<=n; j+=i) visit[j] = true; //标记为非素数 //下面记录素数 int k=0; //统计素数个数 for(int i = 2; i <= n; i++) if(!visit[i]) prime[k++] = i; //存储素数 return k; }
``` - 欧拉筛（线性筛）
-
\(O(n)\)
-
用最小质因子进行筛选，防止对一个数进行重复的处理
int prime[N]; //保存质数，为节约空间，可以适当减小
bool vis[N]; //记录是否被筛
int euler_sieve(int n){ //欧拉筛。返回质数的个数。
int cnt = 0; //记录质数个数
memset(vis,0,sizeof(vis));
memset(prime,0,sizeof(prime));
for(int i=2;i<=n;i++){ //检查每个数，筛去其中的合数
if(!vis[i]) prime[cnt++]=i; //如果没有筛过，是质数，记录。第一个质数是2
for(int j=0; j<cnt; j++){ //用已经得到的质数去筛后面的数
if(i*prime[j] >n) break; //只筛小于等于n的数
vis[i*prime[j]]=1; //关键1。用x的最小质因数筛去x
if(i%prime[j]==0) break; //关键2。如果不是这个数的最小质因子，结束
}
}
return cnt; //返回小于等于n的质数的个数
}
-
对于 \(x\) ，我们遍历到质数表中的 \(p\) ，且发现 \(p|x\) 时，就应当停止遍历质数表。因为：设 \(x=pr(r\ge p)\) （ \(p\) 本身是 \(x\) 的最小质因数），那么对于任意 \(p'>p\) ，有 \(p'x=pp'r=p(p'r)\) ，说明 \(p'x\) 的最小质因数不是 \(p'\) ，我们不应该在此划掉它。
-
使用欧拉筛求解每个元素的最小质因子
-
c++ int prime[N]; //记录质数 int vis[N]; //记录最小质因子 int euler_sieve(int n){ int cnt=0; memset(vis,0,sizeof(vis)); memset(prime,0,sizeof(prime)); for(int i=2;i<=n;i++){ if(!vis[i]){ vis[i]=i; prime[cnt++]=i;} //vis[]记录最小质因子 for(int j=0; j<cnt; j++){ if(i*prime[j] >n) break; vis[i*prime[j]] = prime[j]; //vis[]记录最小质因子 if(i%prime[j]==0) break; } } return cnt; }
质因数分解¶
-
算数基本定理：任何一个数可以拆分为质数的乘积
-
\(n=p_1^{c_1}p_2^{c_2}...\)pi为素数
-
因数个数定理：n的因数个数为\(f(n)=\prod_{i=1}^m(c_i+1)\)
- 可以使用数论分分块计算\(f(n)=\sum_{i=1}^n\lfloor \frac{n}{i} \rfloor\)
-
试除法
-
\(O(\sqrt{n})\)
-
除去所有的质因数后如果剩下的元素大于1，那么这也是一个质因数（最大的）
-
c++ int p[20]; //p[]记录因子，p[1]是最小因子。一个int数的质因子最多有十几个 int c[40]; //c[i]记录第i个因子的个数。一个因子的个数最多有三十几个 int factor(int n){ int m = 0; for(int i = 2; i <= sqrt(n); i++) if(n%i == 0){ p[++m] = i, c[m] = 0; while(n%i == 0) n/=i, c[m]++; //把n中重复的因子去掉 } if(n>1) p[++m] = n, c[m] = 1; //没有被除尽，是素数 return m; //共m个 } -
区间元素质因数分解（由于内存限制，范围不宜太大，1e6）
-
使用欧拉筛，得到每个元素的最小质因数，每个数只需要去除这个最小质因数，然后去找下一个质因数即可
-
预处理复杂度\(O(n)\)，单次查询\(O(logn)\)
-
c++ const int N= 1e6 + 5; int minp[N],prime[N]; void init_prim() { for (int i = 2; i < N; i++) minp[i] = i; int cnt = 0; for (int i = 2; i < N; i++) { if (minp[i] == i) prime[cnt++] = i; for (int j = 0; j < cnt && prime[j] * i < N; j++) { minp[prime[j] * i] = prime[j]; if (i % prime[j] == 0) break; } } } //融合试除法，num较大时也不会出错 unordered_map<int,int> cal_prim(int num){ unordered_map<int,int> res; for(int i=2;i<sqrt(num)&&num>=N;i++){ while(num%i==0){ res[i]++; num/=i; } } if(num>=N){ res[num]++; return res; } while (num != 1) { int p = minp[num]; res[p]++; if(p==0) cout<<"!!!"; num /= p; } return res; } -
pollard_rho启发式法
-
期望复杂度\(O(n^{\frac14}\log n)\)
-
随机算法，猜测、验证因数
-
c++ typedef long long ll; ll Gcd (ll a,ll b){return b? Gcd(b, a%b):a;} ll mult_mod (ll a,ll b,ll n){ //返回(a*b) mod n a %= n, b %= n; ll ret=0; while (b){ if (b&1){ ret += a; if (ret >= n) ret -= n; } a<<=1; if (a>=n) a -= n; b>>=1; } return ret; } ll pollard_rho (ll n){ //返回一个因子，不一定是质因子 ll i=1, k=2; ll c = rand()%(n-1)+1; ll x = rand()%n; ll y = x; while (true){ i++; x = (mult_mod(x,x,n)+c) % n; //(x*x) mod n ll d = Gcd(y>x?y-x:x-y, n); //重要：保证gcd的参数大于等于0 if (d!=1 && d!=n) return d; if (y==x) return n; //已经出现过，直接返回 if (i==k) { y=x; k=k<<1;} } } void findfac (ll n){ //找所有的素因子 if (miller_rabin(n)) { //用miller_rabin判断是否为素数 factor[tol++] = n; //存素因子 return; } ll p = n; while (p>=n) p=pollard_rho(p); //找到一个因子 findfac(p); //继续寻找更小的因子 findfac(n/p); } - 完全平方数的质因数的幂都是偶数
- 两个数相乘为完全平方数要求两个数为奇数的质因子相同，因此对每个数所有奇数的质因子相乘，这个数相同的两个数乘积才会wei
威尔逊定理¶
- 若p为素数则\(p|(p-1)!+1\)，即
- \((p-1)!modp\equiv p-1\)
- 计算\(q!modp\) p是一个素数，q是小于p的最大素数
- \(q!(q+1)...(p-1)\mod{p}=p-1\)
- \(q!\mod{p}=\frac{1}{(q+1)(q+2)...(p-2)}\mod{p}\)
- 分数取模：使用费马小定理\((b/a)\mod{m}=(ba^{m-2})\mod={m}\)
积性函数¶
- 定义在所有正整数上的函数称为算术函数
- 若对于互素的pq有\(f(pq)=f(p)f(q)\)称为积性函数
-
若对任两个数都成立，被称为完全积性函数
-
性质：
- 积性函数的和也是积性函数
欧拉函数¶
- \(\varphi(x)\)表示小于等于x的，与x互素的数的数目
- 当要求连续数的欧拉函数时可以将问题转化为对欧拉函数的计算，实现优化
- 欧拉函数是一个积性函数，对于互素的pq有，\(\varphi(pq)=\varphi(p)\varphi(q)\)
- 性质：\(n=\sum_{d|n}\varphi(d)\)
欧拉函数的计算¶
-
通解公式\(\varphi(n)=n\prod_{i=1}^k(1-\frac{1}{p_i})\)
-
n为素数则\(\varphi(n)=n-1\)
-
\(n=p^k\) p是素数 \(\varphi(n)=p^{k-1}\varphi(p)=p^k-p^{k-1}\)
-
也就是说欧拉函数的求解本质上就是分解质因数，单个欧拉函数的求解为\(O(\sqrt{n})\)
-
c++ int euler(int n){ int ans = n; for(int p = 2; p*p <= n; ++ p){ //试除法：检查从2到sqrt(n)的每个数 if(n%p == 0){ //能整除，p是一个因子，而且是质因子，请思考 ans = ans/p*(p-1); //求欧拉函数的通式 while(n%p == 0) //去掉这个因子的幂，并使得下一个p是质因子 n /= p; //减小了n } } if(n != 1) ans = ans/n*(n-1); //情况(1)：n是一个质数，没有执行上面的分解(或者是分解之后还剩下了一个质数) return ans; } -
欧拉筛计算前n项
-
积性函数都可以使用欧拉筛加速求解
-
\(O(n)\)
-
```c++ #include
using namespace std; const int N = 50000; int vis[N]; //记录是否被筛；或者用于记录最小质因子 int prime[N]; //记录质数 int phi[N]; //记录欧拉函数 int sum[N]; //计算欧拉函数的和 void get_phi(){ //模板：求1～N范围内的欧拉函数 phi[1]=1; int cnt=0; for(int i=2;i<N;i++) { if(!vis[i]) { vis[i]=i; //vis[i]=1; //二选一：前者记录最小质因子，后者记录是否被筛 prime[cnt++]=i; //记录质数 phi[i]=i-1; //情况(1)：i是质数，它欧拉函数值=i-1 } for(int j=0;j
int n; scanf("%d",&n); if(n==1) printf("0\n"); else printf("%d\n",2sum[n-1]+1); return 0; }```
例题¶
- P2568 GCD
- P2398 GCD SUM
- 另一种方法
例题¶
- P2613 【模板】有理数取余
- 也可以使用费马小定理解决
- P2480 古代猪文
- 难，综合
线性代数¶
矩阵¶
-
矩阵乘法：
-
m*n n*u矩阵\(c[i][j]=\sum_{k=1}^{n}A[i][k]*B[k][j]\)
-
c++ for(int i=1;i<=m;i++){ for(int j=1;j<=n;j++){ for(int k=1;k<=1;k++) c[i][j]+=a[i][k]+b[k][j]; } } -
具有结合律分配律，但不满足交换律
矩阵快速幂¶
-
c++ struct matrix{ int m[N][N]; }; //定义矩阵，常数N是矩阵的行数和列数 matrix operator * (const matrix& a, const matrix& b){ //重载*为矩阵乘法。注意const matrix c; memset(c.m, 0, sizeof(c.m)); //清零 for(int i=0; i<N; i++) for(int j=0; j<N; j++) for(int k = 0; k<N; k++) //c.m[i][j] += a.m[i][k] * b.m[k][j]; //不取模 c.m[i][j] = (c.m[i][j] + a.m[i][k] * b.m[k][j]) % mod;//取模 return c; } matrix pow_matrix(matrix a, int n){ //矩阵快速幂，代码和普通快速幂几乎一样 matrix ans; memset(ans.m,0,sizeof(ans.m)); for(int i=0;i<N;i++) ans.m[i][i] = 1; //初始化为单位矩阵，类似普通快速幂的ans=1 while(n) { if(n&1) ans = ans * a; //不能简写为ans *= a，这里的*重载了 a = a * a; n>>=1; } return ans; } -
复杂度\(O(N^3\log_2n)\)，N*N的方阵，求n次幂
加速递推¶
-
加速线性递推问题
-
斐波那契
- \(\begin{pmatrix} F_n & F_{n-1}\end{pmatrix}=\begin{pmatrix} F_{n-1} & F_{n-2}\end{pmatrix}*\begin{pmatrix} 1 & 1\\ 1 & 0\end{pmatrix}\)
- 矩阵幂的和\(A+A^2...\)
- \(\begin{pmatrix} F_n & E\end{pmatrix}=\begin{pmatrix} F_{n-1} & E\end{pmatrix}*\begin{pmatrix} A & 0\\ A & E\end{pmatrix}\)
- P5337 甲苯先生的字符串
- 2851. 字符串转换
路径问题¶
- 解决允许重复经过的指定路径长度的问题
- 使用邻接矩阵M存图，\(G=M^n\)中G[i][j]就表示从i到j经过n条边的路径总数
- 如果是求解最短路径需要对矩阵“幂”进行推广\(M*M=\min_{a=1}^N[M[i][a]+M[a][j]]\)
- P3758 可乐
高斯消元¶
- 用于求解线性方程组
- 参数矩阵在左，右侧添加灯饰右边的数，构成增广矩阵（m*(n+1)）
- 初等变换：
- 交换两行
- 一行乘以非零数k
- 把某一行加到另一行上
- 将左侧矩阵变化为单位矩阵，右侧就是方程的解
- 高斯消元结束后，若存在系数为0，常数不为0的行，则方程无解
- 若系数不为0的行有k个，则说明主元有k个，自由元有n-k个，方程多解
高斯约旦消元法¶
-
过程：
-
从第一列开始选择一个最大的非零系数作为主元，将这一行换到第一行，通过变化2将其转化为1然后消去其他行
-
重复过程，对每一列重复操作
-
在消元过程中设计除法运算，出现浮点数，会造成一定的计算误差
-
复杂度\(O(n^3)\)(三个变量三个方程)
-
```c++ #include
using namespace std; double a[105][105]; double eps = 1e-7; int main(){ int n; scanf("%d",&n); for(int i=1;i<=n;++i) for(int j=1;j<=n+1;++j) scanf("%lf",&a[i][j]); for(int i=1;i<=n;++i){ //枚举列 int max=i; for(int j=i+1;j<=n;++j) //选出该列最大系数，真实目的是选一个非0系数 if(fabs(a[j][i])>fabs(a[max][i])) max=j; for(int j=1;j<=n+1;++j) swap(a[i][j],a[max][j]); //移到前面 if(fabs(a[i][i]) < eps){ //对角线上的主元系数等于0,说明没有唯一解（多解、无解） puts("No Solution"); return 0; } for(int j=n+1;j>=1;j--) a[i][j]= a[i][j]/a[i][i]; //把这一行的主元系数变为1 for(int j=1;j<=n;++j){ //消去主元所在列的其他行的主元 if(j!=i) { double temp=a[j][i]/a[i][i]; for(int k=1;k<=n+1;++k) a[j][k] -= a[i][k]*temp; } } } for(int i=1;i<=n;++i) printf("%.2f\n",a[i][n+1]); //最后得到简化阶梯矩阵 return 0; }
```
-
可以处理多解的模板
-
c++ int Gauss(){//对矩阵进行初步化简，尽可能进行消元，并将自由元放在矩阵的底部，方便处理 int r=0,cnt=0; //cnt表示自由变元个数 for(int c=0;r<equ&&c<var;++r,++c){ int Maxr=r; for(int i=r+1;i<equ;++i) if(abs(a[i][c])>abs(a[Maxr][c])) Maxr=i; if(Maxr!=r){ for(int i=c;i<var+1;++i) swap(a[Maxr][i],a[r][i]); } if(!a[r][c]){ --r;//保持行不变，但列到下一列（即把自由元都放在了末尾） free_xx[cnt++]=c;//存储自由元 continue; } for(int i=r+1;i<equ;++i){ if(!a[i][c]) continue; for(int j=c;j<var+1;++j) a[i][j]^=a[r][j]; } } for(int i=r;i<equ;++i) if(a[i][var]) return -1; //无解 return var-r; //返回自由变元的个数,cnt=var-r }
例题¶
-
在原矩阵右加一个单位矩阵，通过初等变换将左侧变为单位矩阵，右侧就得到逆矩阵
-
做差消除二次项
-
广义线性方程（不再是普通的加法）+ 对多解情况的处理
-
poj1830（方案数目）
- ```c++ const int MAXL(50); const int INF(0x3f3f3f3f); const int mod(1e9+7); int dir[4][2]= {{-1,0},{1,0},{0,1},{0,-1}}; int a[MAXL+50][MAXL+50];//增广矩阵 int x[MAXL+50];//解集 bool free_x[MAXL];//标记是否是不确定的变元
int Gauss(int equ,int var) { memset(x,0); int maxr,col; int k; for(k=0,col=0;k<equ&&col<var;k++,col++) { maxr=k; for(int i=k+1;i
-
POJ1681（求最小方案）
-
自己建立开关间的关系
-
使用异或代替加法，并且存在多解（五穷解），需要通过搜索确定自由元，对自由元的状态枚举，确定最小结果
-
n*n元素方程组
-
-
```c++ #include
#include #include #include #include using namespace std;
const int maxn=250; const int inf=0x3f3f3f3f; int T,n,equ,var,a[maxn][maxn],x[maxn],free_xx[maxn]; int ans; char s[20];
void init(){ //初始化 memset(a,0,sizeof(a)); for(int i=0;i<n;++i){ for(int j=0;j
int Gauss(){//对矩阵进行初步化简，尽可能进行消元，并将自由元放在矩阵的底部，方便处理 int r=0,cnt=0; //cnt表示自由变元个数 for(int c=0;r<equ&&c<var;++r,++c){ int Maxr=r; for(int i=r+1;i
int solve(int t){//尝试遍历自由元方案获取最小解 ans=inf; for(int i=0;i<(1<
int main(){ scanf("%d",&T); while(T--){ scanf("%d",&n); init(); equ=var=n*n; for(int i=0;i<n;++i){ scanf("%s",s); for(int j=0;j<n;++j) if(s[j]'y') a[in+j][nn]=0; else a[in+j][nn]=1; } int t=Gauss(); if(t-1) printf("inf\n"); else printf("%d\n",solve(t)); } return 0; } ```
-
-
hdu5755
- 使用模3加法
线性基¶
- 即向量空间的基向量，可以用于对问题进行压缩
求解¶
-
使用高斯消元法可以求解线性基
-
如果集合中所有元素的二进制位数不同，那么它就是自己的线性基
-
假设最大二进制位数为m，将n个数写成n*m的矩阵，目标就是将其转化为简化行阶梯矩阵（部分是单位矩阵），非零行就是线性基
-
首个非零元素都是1、且是所在列唯一非零元素的行阶梯形矩阵
-
\(\begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix} ， \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{pmatrix} ， \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\\ \end{pmatrix}\)
-
不使用高斯消元法构造(常用)：
-
为了使加入线性基的元素之间的长度各不相同，可以在加入时与已经在线性基的元素中异或实现
-
c++ //不用高斯消元求线性基 #include<bits/stdc++.h> using namespace std; typedef long long ll; const int M=63; ll p[M]; //线性基 bool zero; void Insert(ll x){ for(int i=M;i>=0;i--) if(x>>i == 1) //x的最高位 if(p[i]==0){ p[i]=x; return; } //P[i]还没有，直接让P[i] = x else x^=p[i]; // P[i]已经有了，逐个异或 zero = true; //A有异或和为0的组合 } ll qmax(){ ll ans = 0; for( int i=M;i>=0;i--) ans = max(ans,ans^p[i]); return ans; } int main(){ ll x; int n; scanf("%d",&n); for(int i=1;i<=n;i++) scanf("%lld",&x), Insert(x); printf("%lld\n",qmax()); return 0; } -
高斯消元法(求第k大异或)：
-
```c++ #include
#define N 10100 using namespace std; typedef long long ll; int n; bool zero; //消元后是否产生全0的行 ll a[N]; void Gauss(){ //高斯消元求线性基 int i,k=1; //k标记当前第几行 ll j = (ll)1<<62; //注意不是63 for(;j;j>>=1){ for(i=k;i<=n;i++) if(a[i]&j) break; //找到第j位是1的a[] if(i > n) continue; //没有第j位是1的a[] swap(a[i],a[k]); //把这一行换到上面 for(i=1;i<=n;i++) //生成简化阶梯矩阵(将其他行化为0) if(i != k && a[i]&j) a[i]^=a[k]; k++; } k--; if(k!=n) zero = true; else zero = false; n = k; //线性基中元素的个数 } ll Query(ll k){ //第k小异或和 ll ans=0; if(zero) k--;//存在0 if(!k) return 0; for(int i=n;i;i--){//按二进制位取元素 if(k&1) ans^=a[i]; k >>= 1; } if(k) return -1; return ans; } int main(){ int cnt=0; int T; cin>>T; while(T--){ printf("Case #%d:\n",++cnt); cin>>n; for(int i=1;i<=n;i++) scanf("%lld",&a[i]); Gauss(); int q; cin>>q; while(q--){ ll k; scanf("%lld",&k); printf("%lld\n", Query(k) ); } } } ```
-
注意：
(ll)1<<j如果j>30，即使j是ll，也要将1声明为ll
应用¶
异或空间线性基¶
-
时间复杂度\(O\log{n}\)
-
应用于异或空间，高效求解异或问题
-
假设有n个数，任选组合，使其异或值最大
- 虽然n个数中任选n个数的组合数非常多，但是实际上结果种类就是\(2^m\)(m为最长二进制位数，约为lognmax)，这就从2n缩小到2m
-
可以计算原集合的线性基，对线性基的组合求解
-
性质：
-
在原数组上进行异或计算与在线性基上的计算结果相同
- 是满足以上条件的最小集合
- 线性基中不存在异或为0的组合（这说明有元素可以被其他元素替代）
-
线性基中每个元素的二进制位数均不同
-
最小异或和：
-
线性基中的最小元素
-
最大异或和：
-
如果异或一个基能更大，那就异或（高斯消元法直接异或所有不为0的）
-
第k大异或和：
-
高斯消元法比较方便，按最高位从小到大为数的二进制的每一位，根据二进制取出异或即可（0001-0010-0011...）
-
按照从大到小排序，逐个将元素插入到线性基，能加入就加
-
不适合使用高斯消元法（因为尝试构建梯形，而这题并不重要，只关心能不能加入，加在什么位置不重要）
-
c++ void f() { for (int x = 1; x <= n; x++) { ll t=a[x].n; for (int i = 62; i >= 0; i--) { if ((t >> i) & 1) // 从左向右判断每一位的1 { if (d[i + 1] == 0) // 表示此位的1消不掉 { d[i + 1] = t; res += a[x].v; break; } else t ^= d[i + 1]; // 保存异或结果 } } } }
计算几何¶
基础¶
-
几何坐标通常为实数，使用double类型（%lf）
-
圆周率
const double pi = acos(-1.0); -
符号判断(=0)
-
c++ const int eps=1e-8; int sgn(double x){ if(fabs(x) < eps) return 0; else return x<0?-1:1; } -
大小
-
c++ int dcmp(double x, double y){ //比较两个浮点数 if(fabs(x - y) < eps) return 0; //x==y，返回0 else return x<y ?-1:1; //x<y返回-1，x>y返回1 } -
点的表示
-
c++ struct Point{ double x,y; Point(){} Point(double x,double y):x(x),y(y){} Point operator + (Point B){return Point(x+B.x,y+B.y);} Point operator - (Point B){return Point(x-B.x,y-B.y);} Point operator * (double k){return Point(x*k,y*k);}//与实数相乘 Point operator / (double k){return Point(x/k,y/k);} bool operator == (Point B){return sgn(x-B.x)==0 && sgn(y-B.y)==0;}//近似相等 }; -
两点距离
double Distance(Point A, Point B){ return hypot(A.x-B.x,A.y-B.y); }- 一个计算距离的库函数
-
向量
-
为了简化，用点表示向量，即从(0,0)指向(x,y)的有向线段
-
typedef Point Vector; -
点乘
-
double Dot(Vector A,Vector B){ return A.x*B.x + A.y*B.y; } -
可以使用点乘判断夹角的类型
-
求解长度
c++ double Len(Vector A){return sqrt(Dot(A,A));} double Len2(Vector A){return Dot(A,A);}
-
求夹角
double Angle(Vector A,Vector B){return acos(Dot(A,B)/Len(A)/Len(B));} -
叉乘
-
double Cross(Vector A,Vector B){return A.x*B.y - A.y*B.x;} -
判断向量相对位置
- \(A \times B>0\)：B在A的逆时针方向
- \(A \times B=0\)：共线（方向不确定）
- 检查平行
bool Parallel(Vector A, Vector B){return sgn(Cross(A,B)) == 0;} - \(A \times B<0\)：B在A的顺时针方向
-
计算两向量构成的平行四边形的有向面积(如ABC三点，以A为公共点，AB、AC两边构成平行四边形)
double Area2(Point A,Point B,Point C){ return Cross(B-A, C-A);}
-
向量旋转
-
旋转\(\Theta\)度，有\(x^
=x\cos\Theta-y\sin\Theta \\ y^=x\sin\Theta+y\cos\Theta\) -
c++ Vector Rotate(Vector A, double rad){ return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad)); } -
求法向量
Vector Normal(Vector A){return Vector(-A.y/Len(A), A.x/Len(A));}逆时针旋转90度并取单位向量
-
二维几何¶
点和线¶
-
直线的表示
-
用两个点表示
-
一般式子\(ax+by+c=0\)
-
斜截式\(y=kx+b\)
-
点加斜率
-
点加向量\(P=P_0+vt_0\)，v可以表示\(v=B-A\)，t取不同值可以表示线段、直线、射线
-
c++ struct Line{ Point p1,p2; //（1）线上的两个点 Line(){} Line(Point p1,Point p2):p1(p1),p2(p2){} Line(Point p,double angle){ //（4）根据一个点和倾斜角 angle 确定直线,0<=angle<pi p1 = p; if(sgn(angle – pi/2) == 0){p2 = (p1 + Point(0,1));} else{p2 = (p1 + Point(1,tan(angle)));} } Line(double a,double b,double c){ //（2）ax+by+c=0 if(sgn(a) == 0){ p1 = Point(0,-c/b); p2 = Point(1,-c/b); } else if(sgn(b) == 0){ p1 = Point(-c/a,0); p2 = Point(-c/a,1); } else{ p1 = Point(0,-c/b); p2 = Point(1,(-c-a)/b); } } }; -
点和直线的关系
-
直线上、左、右
-
c++ int Point_line_relation(Point p, Line v){ int c = sgn(Cross(p-v.p1,v.p2-v.p1)); if(c < 0)return 1; //1：p在v的左边 if(c > 0)return 2; //2：p在v的右边 return 0; //0：p在v上 } -
点和线段的关系
-
用叉积判断是否共线，用点积判断是否与端点连线夹角是否180度
-
c++ bool Point_on_seg(Point p, Line v){ //点和线段：0 点不在线段v上；1 点在线段v上 return sgn(Cross(p-v.p1, v.p2-v.p1)) == 0 && sgn(Dot(p – v.p1,p- v.p2)) <= 0; } -
点到直线的距离
-
叉乘求平行四边形面积，除以底边
-
c++ double Dis_point_line(Point p, Line v){ return fabs(Cross(p-v.p1,v.p2-v.p1))/Distance(v.p1,v.p2); } -
点到线段的距离
-
如果点到直线的投影在线段上，那就是最短距离，否则是到两个端点距离的最小值
-
通过点乘法判断与线段成角是否是锐角，判断能否直接取垂直距离
-
c++ double Dis_point_seg(Point p, Segment v){ if(sgn(Dot(p- v.p1,v.p2-v.p1))<0 || sgn(Dot(p- v.p2,v.p1-v.p2))<0) return min(Distance(p,v.p1),Distance(p,v.p2)); return Dis_point_line(p,v); } -
点的投影
-
c++ Point Point_line_proj(Point p, Line v){ double k = Dot(v.p2-v.p1,p-v.p1)/Len2(v.p2-v.p1); return v.p1+(v.p2-v.p1)*k; } -
对称点
-
先求投影点，自然得到对称点
-
c++ Point Point_line_symmetry(Point p, Line v){ Point q = Point_line_proj(p,v); return Point(2*q.x-p.x,2*q.y-p.y); } -
直线的位置关系
-
c++ int Line_relation(Line v1, Line v2){ if(sgn(Cross(v1.p2-v1.p1,v2.p2-v2.p1)) == 0){ if(Point_line_relation(v1.p1,v2)==0) return 1; //1 重合 else return 0; //0 平行 } return 2; //2 相交 } -
直线的交点
-
通过三角形面积比例求出
-
c++ //Line1:ab, Line2:cd Point Cross_point(Point a,Point b,Point c,Point d){ double s1 = Cross(b-a,c-a); double s2 = Cross(b-a,d-a); //叉积有正负 return Point(c.x*s2-d.x*s1,c.y*s2-d.y*s1)/(s2-s1); } -
判断线段是否相交
-
如果一个线段的两端在另一条线段的两端，那么另一个端点与另一线段产生的叉积相反，因此本体两个线段都在对方两端，说明线段相交
-
c++ //Line1:ab, Line2:cd //1相交；0不相交 bool Cross_segment(Point a,Point b,Point c,Point d){ double c1 = Cross(b-a,c-a),c2=Cross(b-a,d-a); double d1 = Cross(d-c,a-c),d2=Cross(d-c,b-c); return sgn(c1)*sgn(c2) < 0 && sgn(d1)*sgn(d2) < 0; } -
判断直线和线段相交
-
c++ bool Cross_segment(Line l, Point a, Point b) { double c1 = Cross(l.p2 - l.p1, a - l.p1), c2 = Cross(l.p2 - l.p1, b - l.p1); return sgn(c1) * sgn(c2) <= 0; } -
线段的交点
-
先判断是否相交，如果相交，使用直线求交点的方法即可。
例题¶
- P1355 神秘大三角 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)
- 点与三条边叉乘（大于零时r++；小于时l++）
- 若l==3||r==3说明在内部
- 若l>0&&r>0说明在外部
- 若l+r==1说明在顶点
- 否则在边上
- P1863 独眼兔
- 叉乘判断方向，点成获取角度，排序
多边形¶
-
点和多边形的关系
-
从P引出射线，转动与每一条边相交，从而判断点在多边形内部还是在多边形外部
-
叉乘计算p在直线的哪一侧，uv判断过p的水平线是否穿过多边形对应的边
-
只使用叉乘判断是否在所有边的同一侧只针对凸多边形有效，本方法对所有多边型都可以使用
-
c++ int Point_in_polygon(Point pt,Point *p,int n){ //点pt，多边形Point *p for(int i = 0;i < n;i++){ //3：点在多边形的顶点上 if(p[i] == pt) return 3; } for(int i = 0;i < n;i++){ //2：点在多边形的边上 Line v=Line(p[i],p[(i+1)%n]); if(Point_on_seg(pt,v)) return 2; } int num = 0; for(int i = 0;i < n;i++){ int j = (i+1)% n; int c = sgn(Cross(pt-p[j],p[i]-p[j])); int u = sgn(p[i].y – pt.y); int v = sgn(p[j].y – pt.y); if(c > 0 && u < 0 && v >=0) num++; if(c < 0 && u >=0 && v < 0) num--; } return num != 0; //1：点在内部; 0：点在外部 } -
多边形的面积
-
相当于从原点连接所有点，做三角形刨分，分别计算面积，由于叉乘的正负性，多边形之外的部分的面积会被抵消，可以求解非凸多边形
-
c++ double Polygon_area(Point *p, int n){ //Point *p表示多边形 double area = 0; for(int i = 0;i < n;i++) area += Cross(p[i],p[(i+1)%n]); return area/2; //面积有正负，返回时不能简单地取绝对值 } -
多边形的重心
-
进行三角形刨分，计算每个三角形的重心，三角形的重心是顶点坐标的平均值，然后对三角形有向面积加权平均，得到结果
-
c++ struct Point{ double x,y; Point(double X=0,double Y=0){x=X,y=Y;} Point operator + (Point B){return Point (x+B.x,y+B.y);} Point operator - (Point B){return Point (x-B.x,y-B.y);} Point operator * (double k){return Point (x*k,y*k);} Point operator / (double k){return Point (x/k,y/k);} }; typedef Point Vector; double Cross(Vector A,Vector B){return A.x*B.y - A.y*B.x;} double Polygon_area(Point *p, int n){ //求多边形面积 double area = 0; for(int i = 0;i < n;i++) area += Cross(p[i],p[(i+1)%n]); return area/2; //面积有正负，不能取绝对值 } Point Polygon_center(Point *p, int n){ //求多边形重心 Point ans(0,0); if(Polygon_area(p,n)==0) return ans; for(int i = 0;i < n;i++) ans = ans+(p[i]+p[(i+1)%n])*Cross(p[i],p[(i+1)%n]); return ans/Polygon_area(p,n)/6; } -
\(/6\)由三部分组成：三角形边中点\(1/2\)，重心在中线\(2/3\)，三角形面积\(1/2\)
-
三角形外心
-
c++ Point circle_center(const Point a, const Point b, const Point c){ Point center; double a1=b.x-a.x, b1=b.y-a.y, c1=(a1*a1+b1*b1)/2; double a2=c.x-a.x, b2=c.y-a.y, c2=(a2*a2+b2*b2)/2; double d =a1*b2-a2*b1; center.x =a.x+(c1*b2-c2*b1)/d; center.y =a.y+(a1*c2-a2*c1)/d; return center; } -
高级知识¶
凸包¶
-
把所有点包含在内的最小凸多边形
-
Andrew算法\(O(nlogn)\)
-
先从最左边的点沿下凸包扫描到最有右边，再从最右边沿上凸包扫描到最左边，合并得到完整凸包。
-
把所有点按照x从小到大进行排序，如果x相同按照y排序，并丢弃相同的点，得到\(p_0,\dots,p_m\)
- 从左向右扫描，求出下凸包，如果一个点在凸包前进方向的左边，则说明在下凸包上，加入到凸包，如果在右边，则说明拐弯了，上个点选择有误，删除最近加入下凸包的点（删除可能发生多次），直至完成。再从右向左得到上凸包。（可以使用叉积判断方侧）
-
-
c++ #include <bits/stdc++.h> using namespace std; const int N = 1e5+1; const double eps = 1e-6; int sgn(double x){ //判断x是否等于0 if(fabs(x) < eps) return 0; else return x<0?-1:1; } struct Point{ double x,y; Point(){} Point(double x, double y):x(x),y(y){} Point operator + (Point B){return Point(x+B.x,y+B.y);} Point operator - (Point B){return Point(x-B.x,y-B.y);} bool operator == (Point B){return sgn(x-B.x) == 0 && sgn(y-B.y) == 0;} bool operator < (Point B){ //用于sort()排序，先按x排序，再按y排序 return sgn(x-B.x)<0 || (sgn(x-B.x)==0 && sgn(y-B.y)<0);} }; typedef Point Vector; double Cross(Vector A,Vector B){return A.x*B.y - A.y*B.x;} //叉积 double Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);} //Convex_hull()求凸包。凸包顶点放在ch中，返回值是凸包的顶点数 int Convex_hull(Point *p,int n,Point *ch){ n = unique(p,p+n)-p; //去除重复点 sort(p,p+n); //对点排序：按x从小到大排序，如果x相同，按y排序 int v=0; //求下凸包。如果p[i]是右拐弯的，这个点不在凸包上，往回退 for(int i=0;i<n;i++){ while(v>1 && sgn(Cross(ch[v-1]-ch[v-2],p[i]-ch[v-1]))<=0) v--; ch[v++]=p[i]; } int j=v; //求上凸包 for(int i=n-2;i>=0;i--){//最右点n-1已经在序列里了，故从n-2开始 while(v>j && sgn(Cross(ch[v-1]-ch[v-2],p[i]-ch[v-1]))<=0) v--; ch[v++]=p[i]; } if(n>1) v--;//最左侧点被重复计入 return v; //返回值v是凸包的顶点数 } Point p[N],ch[N]; //输入点是p[]，计算得到的凸包顶点放在ch[]中 int main(){ int n; cin >> n; for(int i=0;i<n;i++) scanf("%lf%lf",&p[i].x,&p[i].y); int v = Convex_hull(p,n,ch); //返回凸包的顶点数v double ans=0; for(int i=0;i<v;i++) ans += Distance(ch[i],ch[(i+1)%v]); //计算凸包周长 printf("%.2f\n",ans); return 0; }
旋转卡壳¶
- 用两条平行线卡住凸包，可以解决很多问题
- 如果过凸包上的两个点可以画一对平行直线，使凸包上所有点都夹在两条平行线之间或者落在平行线上，那么这两个点称为一对对踵点
- 首先找初始对踵点和平行线（可以取y最大最小的两个），做一条向左，一条向右的水平线，同时逆时针旋转两条线，知道一条线与多边形的一条边重合，这就得到了新的对踵点。重复步骤，知道回到初始的对踵点。
- 如最大距离和就可以通过计算所有对踵点对的距离取最大值得到
最近点对¶
-
给定平面上n个点，求最近的两个点
-
分治法\(O(nlogn)\)
-
按照x坐标排序后划分，直至子集中只有1、2个点，解决子集内的问题
-
求出子集的结果之后进行合并：最近点对在一个子集内部；在不同子集内时，设\(S_1\)中的最短距离为\(d_1\)，\(S_2\)中为\(d_2\)，\(d=min(d_1,d_2)\)，取S1、S2中点划分这样一条区域[mid-d,mid+d]
-
- 将点根据y排序，可以枚举左侧的点，对所有对应的候选点进行计算（实际上最多只有6个）
-
也可以这样考虑，把条形区域内全部点加入并按照y轴排序后，可以枚举一个点，另一个是高[0,dis]的点，这样就能不重复的查找（并且上图保证复杂度，即范围内的点是常数级）
-
c++ const double eps = 1e-8; const int N = 100010; const double INF = 1e20; int sgn(double x){ if(fabs(x) < eps) return 0; else return x<0?-1:1; } struct Point{double x,y;}; double Distance(Point A, Point B){return hypot(A.x-B.x,A.y-B.y);} bool cmpxy(Point A,Point B){ //排序：先对横坐标x排序，再对y排序 return sgn(A.x-B.x)<0 || (sgn(A.x-B.x)==0 && sgn(A.y-B.y)<0); } bool cmpy (Point A,Point B){return sgn(A.y-B.y)<0;} //只对y坐标排序 Point p[N],tmp_p[N]; double Closest_Pair(int left,int right){ double dis = INF; if(left == right) return dis; //只剩1个点 if(left + 1 == right) return Distance(p[left], p[right]);//只剩2个点 int mid = (left+right)/2; //分治 double d1 = Closest_Pair(left,mid); //求s1内的最近点对 double d2 = Closest_Pair(mid+1,right); //求s2内的最近点对 dis = min(d1,d2); int k = 0; for(int i=left;i<=right;i++) //在s1和s2中间附近找可能的最小点对 if(fabs(p[mid].x - p[i].x) <= dis) //按x坐标来找 tmp_p[k++] = p[i]; sort(tmp_p,tmp_p+k,cmpy); //按y坐标排序，用于剪枝。这里不能按x坐标排序 for(int i=0;i<k;i++) for(int j=i+1;j<k;j++){ if(tmp_p[j].y - tmp_p[i].y >= dis) break; //剪枝 dis = min(dis,Distance(tmp_p[i],tmp_p[j])); } return dis; //返回最小距离 } int main(){ int n; cin >>n; for(int i=0;i<n;i++) scanf("%lf%lf",&p[i].x,&p[i].y); sort(p,p+n,cmpxy); //先排序 printf("%.4f\n",Closest_Pair(0,n-1)); //输出最短距离 return 0; }
半平面交¶
-
半平面就是平面的一半，定义一条有向线段的左侧就是代表的半平面，半平面的交一定是凸多边形（可能不闭合）
-
有向直线的定义
-
```c++ struct Line{ Point p; //直线上一个点 Vector v; //方向向量，它的左边是半平面 double ang; //极角，从x正半轴旋转到v的角度 Line(){}; Line(Point p, Vector v):p(p),v(v){ang = atan2(v.y, v.x);} bool operator < (Line &L){return ang < L.ang;} //用于排序 };
```
-
计算半平面交获得的凸多边形
-
\(O(nlogn)\)
- 从逆时针看，所有有向边的斜率是单调递增的，因此先按照斜率单调递增进行排序，然后逐个半平面交，使用双端队列记录，首部为最早加入的半平面，尾部为最新加入的
- 新加入的边所有可能的情况
- 可以直接加入、覆盖队尾、覆盖队首、不能加入
- 覆盖队尾：队尾的两个半平面的交点在新加入边的外面，那么删除队尾半平面
- 覆盖队首：队首的两个半平面的交点在新加入边的外面，那么删除队首半平面
- 不能加入：该情况只能发生在排序的最后一个半平面，如果尾部的半平面交点在首部L1的外面，则要删除
-
注意以上操作都可能在加入一条边时操作多次
-
c++ const double INF = 1e12; const double pi = acos(-1.0); //圆周率，精确到15位小数：3.141592653589793 const double eps = 1e-8; int sgn(double x){ if(fabs(x) < eps) return 0; else return x<0?-1:1; } struct Point{ double x,y; Point(){} Point(double x,double y):x(x),y(y){} Point operator + (Point B){return Point(x+B.x,y+B.y);} Point operator – (Point B){return Point(x-B.x,y-B.y);} Point operator * (double k){return Point(x*k,y*k);} }; typedef Point Vector; double Cross(Vector A,Vector B){return A.x*B.y – A.y*B.x;} //叉积 struct Line{ Point p; Vector v; double ang; Line(){}; Line(Point p,Vector v):p(p),v(v){ang=atan2(v.y,v.x);} bool operator < (Line &L){return ang<L.ang;} //用于极角排序 }; //点p在线L左边，即点p在线L在外面： bool OnLeft(Line L,Point p){return sgn(Cross(L.v,p-L.p))>0;} Point Cross_point(Line a,Line b){ //两直线交点 Vector u=a.p-b.p; double t=Cross(b.v,u)/Cross(a.v,b.v); return a.p+a.v*t; } vector<Point> HPI(vector<Line> L){ //求半平面交，返回凸多边形 int n=L.size(); sort(L.begin(),L.end()); //将所有半平面按照极角排序。 int first,last; //指向双端队列的第一个和最后一个元素 vector<Point> p(n); //两个相邻半平面的交点 vector<Line> q(n); //双端队列 vector<Point> ans; //半平面交形成的凸包 q[first=last=0]=L[0]; for(int i=1;i<n;i++){ //情况1：删除尾部的半平面 while(first<last && !OnLeft(L[i], p[last-1])) last--; //情况2：删除首部的半平面： while(first<last && !OnLeft(L[i], p[first])) first++; q[++last]=L[i]; //将当前的半平面加入双端队列尾部 //极角相同的两个半平面，保留左边：(由于排序只针对斜率，因此可能出现后出现右侧的无用直线) if(fabs(Cross(q[last].v,q[last-1].v)) < eps){ last--; if(OnLeft(q[last],L[i].p)) q[last]=L[i]; } //计算队列尾部半平面交点： if(first<last) p[last-1]=Cross_point(q[last-1],q[last]); } //情况3：删除队列尾部的无用半平面 while(first<last && !OnLeft(q[first],p[last-1])) last--; if(last-first<=1) return ans; //空集 p[last]=Cross_point(q[last],q[first]); //计算队列首尾部的交点。 for(int i=first;i<=last;i++) ans.push_back(p[i]); //复制。 return ans; //返回凸多边形 } int main(){ int T,n; cin>>T; while(T--){ cin>>n; vector<Line> L; L.push_back(Line(Point(0,0),Vector(0,-1))); //加一个半平面F:反向y轴 L.push_back(Line(Point(0,INF),Vector(-1,0))); //加一个半平面E:y极大的向左的直线 while(n--){ double a,b; scanf(“%lf%lf”,&a,&b); L.push_back(Line(Point(0,a),Vector(1,b))); } vector<Point> ans=HPI(L); //得到凸多边形 printf(“%d\n”,ans.size()-2); //去掉人为加的两个点 } return 0; }
圆¶
基本定义¶
-
圆的表示
-
c++ struct Circle{ Point c; //圆心 double r; //半径 Circle(){} Circle(Point c,double r):c(c),r(r){} Circle(double x,double y,double _r){c=Point(x,y);r = _r;} }; -
点和圆的关系
-
c++ int Point_circle_relation(Point p, Circle C){ double dst = Distance(p,C.c); if(sgn(dst – C.r) < 0) return 0; //0 点在圆内 if(sgn(dst – C.r) ==0) return 1; //1 圆上 return 2; //2 圆外 } -
直线和圆的关系
-
根据圆心到直线距离判断
-
c++ int Line_circle_relation(Line v,Circle C){ double dst = Dis_point_line(C.c,v); if(sgn(dst-C.r) < 0) return 0; //0 直线和圆相交 if(sgn(dst-C.r) ==0) return 1; //1 直线和圆相切 return 2; //2 直线在圆外 } -
线段和圆的关系
-
c++ int Seg_circle_relation(Segment v,Circle C){ double dst = Dis_point_seg(C.c,v); if(sgn(dst-C.r) < 0) return 0; //0线段在圆内 if(sgn(dst-C.r) ==0) return 1; //1线段和圆相切 return 2; //2线段在圆外 } -
直线和圆的交点
-
先求圆心在直线上投影，再根据距离和半径求出交点与投影点的距离，结合直线的单位向量得到点的坐标
-
c++ //pa, pb是交点。返回值是交点个数 int Line_cross_circle(Line v,Circle C,Point &pa,Point &pb){ if(Line_circle_relation(v, C)==2) return 0;//无交点 Point q = Point_line_proj(C.c,v); //圆心在直线上的投影点 double d = Dis_point_line(C.c,v); //圆心到直线的距离 double k = sqrt(C.r*C.r-d*d); if(sgn(k) == 0){ //1个交点，直线和圆相切 pa = q; pb = q; return 1; } Point n=(v.p2-v.p1)/ Len(v.p2-v.p1); //单位向量 pa = q + n*k; pb = q - n*k; return 2; //2个交点 }
最小圆覆盖¶
-
找一个最小的圆包围所有的点（点可以在圆上）
-
确定圆有两种情况，由两个点确定（圆心在中点），或者由三个点决定（圆心在外心）
-
算法从一个点开始，每次加入一个新的点，更新最小圆，直到扩展到n个点，设得到的最小圆为 \(C_i\)
-
加入一个点，圆心\(p_1\)，半径\(0\)
-
两个点，圆心 \(p_1p_2\) 中点，半径为线段一半
-
三个点，若 \(p_3\) 在圆上/内，则可以直接忽略这个点，否则 \(p_3\) 一定在圆上，需要重新计算，改为\(p_3\)作为第一个点，然后加入 \(p_1p_2\)（\(p_3\)一定在，从\(p_1p_2\)种找点与\(p_3\)组合）
-
三重循环
-
分别枚举三个点，作为第一二三个点（可能只有两个，不符合才会加入第三个）
-
随机打乱点后期望复杂度\(O(n)\)
-
```c++ #define eps 1e-8 const int N = 1e5+1; int sgn(double x){ if(fabs(x) < eps) return 0; else return x<0?-1:1; } struct Point{ double x, y; }; double Distance(Point A, Point B){return hypot(A.x-B.x,A.y-B.y);} Point circle_center(const Point a, const Point b, const Point c){ Point center; double a1=b.x-a.x, b1=b.y-a.y, c1=(a1a1+b1b1)/2; double a2=c.x-a.x, b2=c.y-a.y, c2=(a2a2+b2b2)/2; double d =a1b2-a2b1; center.x =a.x+(c1b2-c2b1)/d; center.y =a.y+(a1c2-a2c1)/d; return center; } //返回圆心和半径 void min_cover_circle(Point *p, int n, Point &c, double &r){ random_shuffle(p, p + n); //随机函数，打乱所有点。这一步很重要 c=p[0]; r=0; //从第1个点p0开始。圆心为p0，半径为0 for(int i=1;i
int main(){ int n; cin >> n; for(int i=0;i<n;i++) scanf("%lf%lf",&p[i].x,&p[i].y); Point c; double r; //最小覆盖圆的圆心和半径 min_cover_circle(p,n,c,r); printf("%.10f\n%.10f %.10f\n",r,c.x,c.y); return 0; }
```
补充¶
- Pick定理：如果一个多边形的顶点都是格点（整数点），多边形的面积等于边界上格点数目（可以使用线性丢番图求解）的一半加上内部格点数减一\(s=j/2+k-1\)
曼哈顿距离转切比雪夫距离¶
- 由此通过坐标变换，将原先的曼哈顿距离转化为投影的最大长度
-
- 求删除一点后两点距离最大值的最小值
class Solution {
public:
int minimumDistance(vector<vector<int>> &points) {
multiset<int> xs, ys;
for (auto &p : points) {
xs.insert(p[0] + p[1]);
ys.insert(p[1] - p[0]);
}
int ans = INT_MAX;
for (auto &p : points) {
int x = p[0] + p[1], y = p[1] - p[0];
xs.erase(xs.find(x));
ys.erase(ys.find(y));
ans = min(ans, max(*xs.rbegin() - *xs.begin(), *ys.rbegin() - *ys.begin()));
xs.insert(x);
ys.insert(y);
}
return ans;
}
};
进制¶
-
难点，找出n,m使得\(2^m-2^n=k(k>0)\)
- \(2^n(2^{m-n}-1)\)即可以看作一个二的幂减一再左移（形似11110000）
-
c++ int d = abs(a[i] - sum); int p = lowbit(d); int e = d + p; if (__builtin_popcount(e) == 1) { if (a[i] > sum) bit[__lg(e)]++, bit[__lg(p)]--; else bit[__lg(e)]--, bit[__lg(p)]++; } else { cout << "No" << '\n'; return; } -
拆分，统计唯一的给与与接受方案（bit均衡，全零）
-
如果一个点与s的差是2的幂则有两种分割方式
- 先正常统计只有一种分割方式的元素，以及有两种分割方式的元素（记录元素差值的对数，但不直接加入bit）
-
从大到小遍历，由于bit[i]会受到pow[i-1]分割方式的影响，可以通过i均衡确定i-1的分割（不能用i的分割方式来均衡i，否则会影响到i+1无法均衡）
-
for (int i = 30; i >= 0; i--) { bit[i] += (pow1[i] - pow2[i]); if (i == 0) break; if (bit[i] < 0) { pow1[i - 1] -= -bit[i]; bit[i - 1] -= -bit[i]; if (pow1[i - 1] < 0) {//shu'mu cout << "No" << '\n'; return; } } else { pow2[i - 1] -= bit[i]; bit[i - 1] += bit[i]; if (pow2[i - 1] < 0) { cout << "No" << '\n'; return; } } }
补充¶
斐波那契¶
- \(gcd(F(n),F(n-1))=1\)
-
\(F[m+n]=F[m-1]F[n]+F[m]F[n+1]\)
-
\(gcd(F(n),F(m))=F(gcd(n,m))\)
- P1306 斐波那契公约数
斐波那契编码¶
- 将斐波那契数作为“进制”
-
整除¶
- 用d(n)表示可以整除n的元素数目
- 若\(n=p^{a_1}p^{a_2}\dots p^{a_k}\)，则显然有\(d(n)=(a_1+1)(a_2+1)\dots(a_k+1)\)
- 能整除n的元素实际上就是使用n的质因数可以组合出多少个不同的数，而每个数字种每个质因数的数目可以是\(0,1,\dots,a_i\)，共\(a_i+1\)
- 若\(\gcd(a,n)=1\)说明没有相同的素因数，则显然有\(d(an)=d(a)d(n)\)，为积性函数
- Problem - F - Codeforces
例题¶
-
lcm拆分|容斥原理
-
质因数的个数
-
求正整数N(N>1)的质因数的个数。 相同的质因数需要重复计算。如120=2*2*2*3*5，共有5个质因数。//如果不是重复计数，那就是求质因数的种类
#include<iostream>
#include<cmath>
using namespace std;
int main(){
int n;
int ans=0;
while(cin>>n){
ans=0;
for(int i=2;i<=sqrt(n);i++){
while(n%i==0){
n=n/i;//事实上这里的i一定是质数，如果不是质数那么一定可以分解成更小的质数，在这之前就已经被除去了，因此这里除的一定是质数
ans++;
}
}
if(n>1){
ans++;//比n的平方根大的质因数最多只能有一个（n>1就说明有）
}
cout<<ans<<endl;
}
}
-
由于浮点数会损失精度，因此要寻找其他方法表示斜率
-
数学找规律
-
分类讨论
-
进制、信息论
-
- 特殊题目：将循环小数转化为分数
- 对于纯循环小数（从第一位开始循环）
- 设循环节长度为 2 ，设为 a ，即 n=0. aaaaa...
- 那么有 100n=a+n
- 即 \(n=\frac{x}{y}=\frac{a}{99}\) 得到结果
- 即 \(n=\frac{a}{10^{\lfloor lg(a) \rfloor+1}-1}\)
- 混循环小数（循环节前面存在非循环片段）
- 设 n=baaa...
- \(n=\frac{b.aaa\dots}{10^{\lfloor lg(b) \rfloor+1}}=\frac{b(10^{\lfloor lg(a) \rfloor+1}-1)+a}{10^{\lfloor lg(b) \rfloor+1}*(10^{\lfloor lg(a) \rfloor+1}-1)}\)
- 最后再通过 gcd 对 xy 化简
- 对式子重新组织，分离变量，实现降维
- 求解阶乘的逆元
fac[0]=1;
for(int i=1;i<=n;i++) fac[i]=fac[i-1]*i%mod;
invfac[n]=expow(fac[n],mod-2);
//根据最后一项求出每一项的逆元
for(int i=n-1;i>=0;i--) invfac[i]=invfac[i+1]*(i+1)%mod;