MIT-18.06

Posted on 2023-01-26 Edited on 2023-01-30 In 笔记

Lec 4 A = LU

\[ (AB)^{-1} = B^{-1}A^{-1} \\ (A^{T})^{-1} = (A^{-1})^{T}\\ \]

GAMES-202

Posted on 2023-01-06 Edited on 2023-01-18 In 笔记

导论

Graphics (Hardware) Pipeline

OpenGL

A. Place objects/models

B. Set up an easel

C. Attach a canvas to the ease

D. Paint to the canvas

MIT 18.01 单变量微积分

Posted on 2023-01-05 Edited on 2023-02-01 In 笔记

常用公式

三角函数

\[ \sin^2x + \cos ^ 2 x = 1\\ \]

求导

三角函数求导

\[ \sin'x = \cos x\\ \cos'x = - \sin x\\ \tan'x = \frac{1}{\cos^2x} = \sec^2x\\ \arctan' x = \frac{1}{1+x^2}\\ \]

指数函数和对数函数

\[ \frac{\text d}{\text d x}\ln x = \frac{1}{x}\\ \frac{\text d}{\text d x}e^x = e^x\\ \frac{\text d}{\text dx}a^x = \ln a \cdot a^x\\ \frac{\text d}{\text d x} \ln u = \frac{1}{u}\frac{\text d u}{\text d x} \]

导数（derivative)

Lecture1: 导数的几何解释

如图, $P$点的导数便是过$P$点的切线(tangent line)的斜率, 记作$f'(P)$.

什么样子的线是切线? 我们先看由点$PQ$构成的割线(secant line). 当$PQ$两点的距离趋于$0$时, 割线与切线便无限接近.

如图. 可得 \[ f'(x) = \lim_{\Delta x \to 0} \frac{f(x_0 + \Delta x) - f(x_0)}{\Delta x} \]

一些符号

导数(牛顿记法): $f'$

导数(莱布尼兹记法): $\frac{\text df}{\text dx},\frac{\text dy}{\text d x},\frac{\text d}{\text dx}f,\frac{\text d}{\text dx}y$

一个有用的公式(推导见手写笔记): \[ \frac{\text d}{\text d x}x^n = n \cdot x^{n-1} \]

Lecture2: 极限, 连续性和三角函数的极限

我们用$\frac{\Delta y}{\Delta x}$来计量某个函数在一段时间内的平均变化率. 当$\Delta x \to 0$时, 这便变成了瞬时变化率.

Limits and Continuity

简单的极限可以直接带入运算:

但是对于 \[ \lim_{\Delta x \to 0}\frac{\Delta f}{\Delta x} \] 它不是一个简单极限. 因为我们不能将$\Delta x = 0$代入运算, 这样永远会得到$\frac{0}{0}$. 注意: \[ \lim_{x \to x_0} \] 暗示了$x \neq x_0$.

左右极限(Left and Right Limt)

我们用以下两个记号来表示左极限和右极限: \[ \lim_{x\to x_0^-} (左极限, 此时x < x_0)\\ \lim_{x\to x_0^+} (右极限, 此时x > x_0)\\ \]

连续性(Continuity)

我们称$f(x)$在$x_0$连续, 当且仅当 \[ \lim_{x \to x_0}f(x) = f(x_0) \] 这意味着:

$f(x)$在点$x_0$的极限必须存在, 且左右极限相等
$f(x_0)$是有定义的
$\lim_{x\to x_0}f(x)$与$f(x_0)$相等

如图, $\lim_{x \to 0^+}f(x) = 1$, 但是$f(0) = 0$. 因此该函数在点$0$不连续.

以下是不连续的几种情况:

Removable Discontinuity

该点的左右极限存在且相等, 但与该点的函数值不等(或该点的函数值未定义)

Jump Discontinuity

左右极限均存在, 但不相等.

Infinite Discontinuity

$x = 0$时, 左极限为负无穷, 右极限为正无穷.

Other (ugly) discontinuities

导数的图像(Picturing the derivative)

注意导数的图像和原函数的图像没有什么相关性. 而且对一个奇函数求导, 它的导数一定是一个偶函数.

定理: 可导必连续(Theorem: Differentiable Implies Continuous)

\[ 如果f在x_0处可导, 那么f必定在x_0处连续. \]

Lecture 3: 导数的四则运算及三角函数

求导公式(Derivative Formulas)

特殊求导公式

\[ \lim_{\theta \to 0}\frac{\sin\theta}{\theta} = 1\\ \lim_{\theta\to0}\frac{1 - cos\theta}{\theta} = 0 \]

证明:

sin, cos

\[ \frac{\text d}{\text d x}\sin x = \cos x\\ \frac{\text d}{\text d x}\cos x = -\sin x\\ \]

推导过程:

一般求导公式

加法法则 \[ (u+v)' = u' + v' \]
乘积法则

\[ (uv)' = u'v + uv' \]

除法法则(quotient rule) \[ (\frac{u}{v})' = \frac{u'v - uv'}{v^2}(v \neq 0) \]
证明: 加法法则
证明: 乘法法则

几何上的理解: $(uv)' \cdot \Delta x = (u + \Delta u)(v + \Delta v) - uv$, 即图中粉色, 白色, 黄色部分之和. 又因为$\Delta u \cdot \Delta v \approx 0$, 所以有 \[ (uv)' \approx \frac{u \Delta v + v \Delta u}{\Delta x} \\ = uv’ +u'v \]

证明: 除法法则

Lecture 4:链式法则和高阶导数

链式法则(Chain Rule)

\[ \frac{\text{d}y}{\text{d}t} = \frac{\text{d} y}{\text{d}x} \cdot \frac{\text{d}x}{\text{d}t} \]

高阶导数(Higher Derivatives)

Lecture 5 隐函数微分及逆函数导数

隐函数微分(Implicit Differentiation)

有些函数(例如$x^2 + y^2 = 1$)不具有明显的$y = f(x)$形式, 但是能通过变形得到(本例可得$y = \sqrt{1 - x^2}$). 这样的函数称为隐函数. 对隐函数求导, 可在方程两边同时对$x$求导, 得到一个包含$\frac{\text{d}y}{\text{d}x}$的式子, 再求解$\frac{\text{d}y}{\text{d}x}$即可.

反函数(Inverse Function)

若$y = f(x),g(x) = y$, 则$g$是$f$的反函数, 记作$f^{-1}$. 函数和它的反函数关于$y = x$对称.

Lecture 6 指数函数与对数函数，对数微分，双曲线函数

指数函数，对数函数和自然对数

见课程pdf和手写笔记

### 对数微分(Logarithmic Differentiation)

Lecture 7 连续性和第一次考试复习

导数伪装的极限

(普林斯顿微积分读本修订版，人民邮电出版社，ISBN: 978-7-115-43559-0, P101)

Leture 9 线性和二阶近似

线性近似(Linear Approximation)

\[ f(x) \approx f(x_0) + f'(x_0)(x - x_0) \]

注意上表近似符号左边和右边的函数. 如果$f(x_0) + f'(x)(x - x_0)$比$f(x)$更容易计算, 便可以用线性近似来加速计算.

二阶近似(Quadratic Approximations)

一阶(线性)近似可能精度不够, 这时便需要二阶或更高阶的近似. \[ f(x) \approx f(x_0) + f'(x_0)(x - x_0) + \frac{f''(x_0)}{2}(x-x_0)^2 \]

关于为什么二阶导的系数是$\frac{1}{2}$:

可以看出, 二阶近似的效果更好:

基础的结论:

Lecture 10: 曲线构图

Lecture 10: Curve Sketching

利用$f'$和$f''$的符号来画出$f$的近似图像

$f' > 0$, $f$递增. 反之亦然
$f'' > 0$, $f'$递增, 反之亦然

画图的一般步骤:

当$f'(x_0) = 0$时, 点$x_0$称为驻点(临界点, critical points).

当$f''(x_0) = 0$时, 点$x_0$称为拐点(inflection point).

[模板]数学

Posted on 2022-08-23 Edited on 2022-09-10 In 模板

[TOC]

答案不取模，亲人两行泪

212370440130137957

(↑这是一个很大的质数)

19260817

(↑这是一颗很大的子弹)

判断同余时用$(a - b) \mod m== 0$

cout不使用科学计数法输出，设置精度:

1 2	cout.setf(ios::fixed,ios::floatfield); cout.precision(2);

# 杂项

## 神秘结论

\[ \sum_{i = 1}^n (d(i))^2 \ll n \log^3 n (d为因子个数) \]

## 大整数

(我也不知道从哪copy)的

struct bign{
    int d[maxn], len;
 
	void clean() { while(len > 1 && !d[len-1]) len--; }
 
    bign() 			{ memset(d, 0, sizeof(d)); len = 1; }
    bign(int num) 	{ *this = num; } 
	bign(char* num) { *this = num; }
    bign operator = (const char* num){
        memset(d, 0, sizeof(d)); len = strlen(num);
        for(int i = 0; i < len; i++) d[i] = num[len-1-i] - '0';
        clean();
		return *this;
    }
    bign operator = (int num){
        char s[20]; sprintf(s, "%d", num);
        *this = s;
		return *this;
    }
 
    bign operator + (const bign& b){
        bign c = *this; int i;
        for (i = 0; i < b.len; i++){
        	c.d[i] += b.d[i];
        	if (c.d[i] > 9) c.d[i]%=10, c.d[i+1]++;
		}
		while (c.d[i] > 9) c.d[i++]%=10, c.d[i]++;
		c.len = max(len, b.len);
		if (c.d[i] && c.len <= i) c.len = i+1;
        return c;
    }
    bign operator - (const bign& b){
        bign c = *this; int i;
        for (i = 0; i < b.len; i++){
        	c.d[i] -= b.d[i];
        	if (c.d[i] < 0) c.d[i]+=10, c.d[i+1]--;
		}
		while (c.d[i] < 0) c.d[i++]+=10, c.d[i]--;
		c.clean();
		return c;
    }
    bign operator * (const bign& b)const{
        int i, j; bign c; c.len = len + b.len; 
        for(j = 0; j < b.len; j++) for(i = 0; i < len; i++) 
			c.d[i+j] += d[i] * b.d[j];
        for(i = 0; i < c.len-1; i++)
            c.d[i+1] += c.d[i]/10, c.d[i] %= 10;
        c.clean();
		return c;
    }
    bign operator / (const bign& b){
    	int i, j;
		bign c = *this, a = 0;
    	for (i = len - 1; i >= 0; i--)
    	{
    		a = a*10 + d[i];
    		for (j = 0; j < 10; j++) if (a < b*(j+1)) break;
    		c.d[i] = j;
    		a = a - b*j;
    	}
    	c.clean();
    	return c;
    }
    bign operator % (const bign& b){
    	int i, j;
		bign a = 0;
    	for (i = len - 1; i >= 0; i--)
    	{
    		a = a*10 + d[i];
    		for (j = 0; j < 10; j++) if (a < b*(j+1)) break;
    		a = a - b*j;
    	}
    	return a;
    }
	bign operator += (const bign& b){
        *this = *this + b;
        return *this;
    }
 
    bool operator <(const bign& b) const{
        if(len != b.len) return len < b.len;
        for(int i = len-1; i >= 0; i--)
            if(d[i] != b.d[i]) return d[i] < b.d[i];
        return false;
    }
    bool operator >(const bign& b) const{return b < *this;}
    bool operator<=(const bign& b) const{return !(b < *this);}
    bool operator>=(const bign& b) const{return !(*this < b);}
    bool operator!=(const bign& b) const{return b < *this || *this < b;}
    bool operator==(const bign& b) const{return !(b < *this) && !(b > *this);}
 
    string str() const{
        char s[maxn]={};
        for(int i = 0; i < len; i++) s[len-1-i] = d[i]+'0';
        return s;
    }
};
 
istream& operator >> (istream& in, bign& x)
{
    string s;
    in >> s;
    x = s.c_str();
    return in;
}
 
ostream& operator << (ostream& out, const bign& x)
{
    out << x.str();
    return out;
}

使用例:

bign GCD(bign a,bign b)
{
	// cout<<"@@@ "<<a<<endl<<"### "<<b<<endl;;
	bign k=0,c;
	while(b>k)
	{
		c=a%b;
		a=b;
		b=c;
	}
	return a;
}
int T;
int vis[110];
vector <bign> prime;
void GetPrime()
{
	for(long long i=2;;i+=1)
	{
		int flag=1;
		for(long long j=2;j<i;j+=1)
			if(i%j==0) flag=0;
		if(flag)
		{
			bign k=i;
			prime.push_back(k);
		}
		if((int)prime.size()>100) break;
	}
}
int main()
{
	cin>>T;
	GetPrime();
	while(T--)
	{
		bign n;
		cin>>n;
		bign a=1,b=1,k=1,sum=1;
		int pos=-1;
		for(int i=0;i<100;i++)
		{
			if(sum*prime[i]<=n) sum=sum*prime[i],pos=i;
			else break;
		}
 
		for(int i=0;i<=pos;i++)
		{
			
			b=b*prime[i];
			a=a*(prime[i]+k);
			
			// cout<<"@@@ "<<s<<endl;
			
		}
		bign s=GCD(a,b);
		b=b/s,a=a/s;
		cout<<b<<"/"<<a<<endl;
	}
	return 0;
 
}

## 随机数发生器

注意不要用Windows下的$rand()$

1
2
3

mt19937 Rand(seed);
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());//cf最好用这个防止被hack
Rand();//返回一个随机数

### 均匀分布

// 包含random头文件
#include<random>
#include<iostream>
int main() {
	//定义均匀分布对象，均匀分布区间（a,b)为(-10，10)
	std::uniform_int_distribution<int> uid{ -10,10 };
	//均匀分布区间(a,b)
	std::cout << uid.a() <<" "<< uid.b() << std::endl;
	//定义随机数种子
	std::random_device rd;
	//定义默认随机数生成器
	std::default_random_engine dre{ rd() };
	std::cout << "五个随机数如下" << std::endl;
	for (int i = 0; i < 5; ++i)
		std::cout << uid(dre) << std::endl;
	return 0;

}

## 计时

clock_t start = clock();
// do something...
clock_t end   = clock();
cout << "花费了" << (double)(end - start) / CLOCKS_PER_SEC << "秒" << endl;

# 数论

## 同余式

同余属于等价关系, 满足自反、对称、传递三大关系 \[ \begin{align} a \equiv a \pmod m\\ a \equiv b \pmod m \iff b \equiv a \pmod m\\ a \equiv b \pmod m \and b \equiv c \pmod m \Rightarrow a \equiv c \pmod m \end{align} \] 加、减、乘与一般运算一致. 若$a \equiv b \pmod m$且$c \equiv d \pmod m$, 则有 \[ a+c \equiv b+d \pmod m\\ a-c \equiv b-d \pmod m\\ a*c \equiv b*d \pmod m \] 除法比较特殊. 若$ac \equiv bc \pmod m$, 则 \[ a \equiv b \pmod{\frac{m}{\gcd(c,m)}} \]

其它性质:

\[ \begin{align} ax \equiv 0 \pmod p \iff x \equiv 0 \pmod{\frac{p}{\gcd(p,a)}} \end{align}\\ x \% b == 0 \and \frac{x}{b} \equiv 0 \pmod p \iff x \equiv 0 \pmod{bm}\\ a \equiv b \pmod m \iff a^n \equiv b^n \pmod m\\ a \equiv b \pmod m \and n|m \Rightarrow a \equiv b \pmod n \]

## 快速幂

计算$x^d \mod p$

注意$a^b % p = ((a % p)^b) % p$

ll Pow(ll x,ll d){
	ll tans = 1;
	if(d == 0)return 1%p;
    x %= p;
	ll a = Pow(x,d/2);
	tans = a*a%p;
	if(d%2)tans = tans*x%p; 
	return tans%p;
}

LL gcd(LL a,LL b)
{
    while(b)
    {
        LL tmp=b;
        b=a%b;
        a=tmp;
    }
    return a;
}

## 一次不定方程

$\sum_{i = 1}^n a_i \cdot x_i = c$有整数解当且仅当$\gcd(a_1,a_2,...,a_n) | c$

## GCD & exGCD

x,y,d为全局变量.求x和y使得 $ax + by = d$ 且 $|x| + |y|$最小.

其中 $d = gcd(a,b)$

ll gcd(ll a,ll b){
	return b?gcd(b,a%b):a;
}
void gcd(ll a,ll b,ll & d,ll & x,ll & y){
	//notice d = gcd(a,b)
	if(!b)d = a,x = 1,y = 0;
	else{
		gcd(b,a%b,d,y,x);
		y -= x*(a/b);
	}
}

## 逆元

返回模p下a的逆.不存在则返回-1.

ll inv(ll a,ll p){
	ll d,x,y;
	gcd(a,p,d,x,y);
	return d == 1?(x+p)%p:-1;
}

### 线性筛逆元:

返回模p下1...n的逆.保存在invs里

void get_inv(){
	inv[1] = 1;
	for(int i = 2;i < maxn;i++){
		inv[i] = (ll)(p - p/i)*inv[p%i]%p;
	}
}

## a,b不互质时

\[ \frac{a}{b} \mod m = \frac{a \mod (mb)}{b} \]

证: \[ \begin{align} &\frac{a}{b} = km + x (x < m)\\ \Rightarrow \space &a =kbm + bx\\ \Rightarrow \space &a \mod (bm) = bx\\ \Rightarrow \space &\frac{a \mod (bm)}{b} = x\\ \Rightarrow \space &\frac{a \mod (bm)}{b} = \frac{a}{b} \mod m \end{align} \]

## 中国剩余定理(CRT)

求同余方程组 \[ \begin{cases} x &\equiv a_1 \pmod {m_1} \\ x &\equiv a_2 \pmod {m_2} \\ &\vdots \\ x &\equiv a_k \pmod {m_k} \\ \end{cases} \] 在$[0,\prod_{i = 1}^nm_i])$中的解. 要求任意$m_i$两两互质

输入: int n 方程组个数, ll a[] ll m[],方程组的参数和模数

复杂度:$O(nlogM)$, 其中$M$与$m_i$同阶

输出: 方程组在$[0,\prod_{i = 1}^nm_i])$中的解.

前置: 逆元

调用 CRT(n,a,m)

ll CRT(int n,ll * a,ll * m){
	ll M = 1;
	for(int i = 1;i <= n;i++)M *= m[i];
	ll ans = 0;
	for(int i = 1;i <= n;i++){
		ll mt = M/m[i];
		ll ret = inv(mt,m[i]);
		ans = (ans + a[i]*mt%M*ret%M)%M;
	}
	return (ans + M)%M;
}

## Eratosthenes筛

void prime(){
	for(int i = 2;i <= n;i++){
		if(!n_prime[i]){
			for(int j = i*2;j <= n;j += i)n_prime[j] = 1;
		}
	}
	n_prime[1] = 1;
}

## 线性筛

筛素因子时记得直接去掉最后一个素因子, 否则复杂度不对

### 线性筛素数

bool not_p[maxn];
vector<int>prime;
void euler_sieve(int n){
	not_p[1] = 1;
	for(int i = 2;i < n;i++){
		if(!not_p[i]){
			prime.push_back(i);
		}
		for(int p:prime){
			if((ll)p*i >= n)break;
			not_p[p*i] = 1;
			if(i%p == 0)break;
		}
	}
}

### 线性筛$\varphi$

$\varphi(n)$: 小于等于$n$的正整数中与$n$互质的数个数

$\varphi$的一些性质(以下$p$代表素数) \[ \varphi(p) = p-1\\ 如果i \% p == 0, \varphi(i*p) = \varphi(i)*p\\ 否则\varphi(i*p) = \varphi(i)*(p-1) \]

phi[n]为$\varphi(n)$,prime里存的素数

int phi[maxn];
bool not_p[maxn];
vector<int>prime;
void eular_sieve(int n){
	not_p[1] = phi[1] = 1;
	for(int i = 2;i < n;i++){
		if(!not_p[i]){
			phi[i] = i-1;
			prime.push_back(i);
		}
		for(int p:prime){
			if((ll)p*i >= n)break;
			not_p[p*i] = 1;
			if(i%p)phi[i*p] = phi[i]*(p-1);
			else{
				phi[i*p] = phi[i]*p;
				break;
			}
		}
	}
}

### 线性筛$\mu$

\[ \mu(n)= \begin{cases} 1&n=1\\ 0&n\text{ 含有平方因子}\\ (-1)^k&k\text{ 为 }n\text{ 的本质不同质因子个数}\\ \end{cases} \]

int phi[maxn],mu[maxn];
bool not_p[maxn];
vector<int>prime;
void sieve_mu(){
	not_p[1] = mu[1] = 1;
	for(int i = 2;i < maxn;i++){
		if(!not_p[i]){
			mu[i] = -1;
			prime.push_back(i);
		}
		for(int p:prime){
			if(ll(p)*i >= maxn)break;
			not_p[p*i] = 1;
			if(i%p)mu[i*p] = -mu[i];
			else{
				mu[i*p] = 0;
				break;
			}
		}
	}
}

## 杜教筛

详细请见杜教筛笔记.

杜教筛快速筛$\varphi$和$\mu$的前缀和.$maxn = n^{\frac{2}{3}}$时复杂度为$O(n^{\frac{2}{3}})$

typedef long long ll;
const int maxn = 1664510;
bool not_p[maxn];
vector<int>prime;
ll mu[maxn],phi[maxn];
void sieve(){
	phi[1] = mu[1] = 1;
	for(int i = 2;i < maxn;i++){
		if(!not_p[i]){
			mu[i] = -1;
			phi[i] = i-1;
			prime.push_back(i);
		}
		for(int p:prime){
			if(1LL*p*i >= maxn)break;
			not_p[i*p] = 1;
			if(i%p){
				phi[i*p] = phi[i]*(p-1);
				mu[i*p] = -mu[i];
			}
			else{
				phi[i*p] = phi[i]*p;
				mu[i*p] = 0;
				break;
			}
 		}
	}
	for(int i = 1;i < maxn;i++)mu[i] += mu[i-1],phi[i] += phi[i-1];
}
map<int,ll>Mu,Phi;
ll du_mu(int x){
	if(x < maxn)return mu[x];
	else if(Mu[x])return Mu[x];
	ll tans = 1;
	for(int l = 2,r = 0;l <= x;l = r+1){
		r = x/(x/l);
		tans -= 1LL*(r - l + 1)*du_mu(x/l);
	}
	return Mu[x] = tans;
}
ll du_phi(int x){
	if(x < maxn)return phi[x];
	else if(Phi[x])return Phi[x];
	ll tans = 1LL*x*(x+1)/2;
	for(int l = 2,r = 0;l <= x;l = r+1){
		r = x/(x/l);
		tans -= 1LL*(r-l+1)*du_phi(x/l);
	}
	return Phi[x] = tans;
}

## 欧拉函数$φ(x)$

ll cal_phi(ll x){
	int m = (ll)sqrt(x+0.5);
	ll tans = x;
	for(int i = 2;i <= m;i++)if(!(x%i)){
		tans = tans/i*(i-1);
		while(!(x%i))x /= i;
	}
	if(x > 1)tans = tans / x * (x-1);
	return tans;
}

### 线性筛$φ(x)$

void sieve_phi(int n){
	phi[1] = 1;
	for(int i = 2;i <= n;i++)if(!phi[i])
		for(int j = i;j <= n;j += i){
			if(!phi[j])phi[j] = j;
			phi[j] = phi[j]/i * (i-1);
		}
}

## 有理数取模

求$a/b$ % $p$ （p为质数)

由费马小定理$b^{p-1} ≡ 1 \pmod p$ 可得 $b^{p-2} ≡ b^{-1} \pmod p$

故只需求$a*b^{p-2} \bmod p$

a,b过大则先膜一下

## 欧拉定理

若$a,p$互质,则有

$a^{φ(p)} ≡ 1 \pmod p$

### 扩展欧拉定理(无需a,p互质)

$b ≥ φ(p)$时:

$a^b ≡ a^{b \space \bmod \space φ(p) +φ(p)} \pmod p$

例题: 给出$w_1,w_2,..w_n$和$q$组询问$l \space r$. 求$w_l^{w_{l+1}^{...^{w_r}}} \pmod p$的值

//CF906D 
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll maxn = 1e5+10;
ll raw[maxn],n,mod;
long long r1aw[maxn];
map<ll,ll>vis;
ll phi(ll x){
	ll m = (ll)sqrt(x+0.5);
	ll tans = x;
	for(ll i = 2;i <= m;i++)if(!(x%i)){
		tans = tans/i*(i-1);
		while(!(x%i))x /= i;
	}
	if(x > 1)tans = tans / x * (x-1);
	return tans;
}
ll Pow(ll x,ll d,ll p){
	ll tans = 1;
	if(d == 0)return 1%p;
    x %= p;
	ll a = Pow(x,d/2,p);
	tans = a*a%p;
	if(d%2)tans = tans*x%p; 
	return tans%p;
}
ll w[maxn],p0[maxn];
int cal(int l,int r,int k) {
	if (l==r) return min(w[l],w[l]%p0[k]+p0[k]);
	if (p0[k]==1) {
		return 1;
	} else {
		int c=cal(l+1,r,k+1);
		if (w[l]==1) return 1;
		ll v=Pow(w[l],c,p0[k]);
		ll d=1ll<<60,e=1;
		for(int i = 0;i < c;i++) { d/=w[l]; e*=w[l]; if (d==0) break; }
		if (d==0) return v+p0[k]; else return min(e,e%p0[k]+p0[k]);
	}
}
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	int n;
	cin >> n >> mod;
	int idx = 0;
	p0[idx] = mod;
	while(p0[idx] != 1){p0[idx+1] = phi(p0[idx]);idx++;}
	for(ll i = 1;i <= n;i++){
		cin >> w[i];
	}
	int q;
	cin >> q;
	while(q--){
		int l,r;
		cin >> l >> r;
		cout << cal(l,r,0)%mod << endl;
	}
	return 0;
}

## 原根

proot(p) 返回素数p的原根

vector<int>pri;
bool g_test(ll g,ll p){
	for(ll x:pri){
		if(Pow(g,(p-1)/x,p) == 1)return 0;
	}
	return 1;
}
ll proot(ll p){
	ll tmp = p-1;
	for(ll i = 2;i <= tmp/i;i++){
		if(tmp%i == 0){
			pri.push_back(i);
			while(tmp%i == 0)tmp /= i;
		}
	}
	if(tmp != 1)pri.push_back(tmp);
	ll g = 1;
	while(1){
		if(g_test(g,p))return g;
		g++;
	}
}

## 离散对数(BSGS)

### 原版

只支持模数为$p$的情况.

bsgs(ll a,ll b,ll p)返回最小的$x$使得$a^x \equiv b\pmod{p}$, 不存在则返回$-1$. 要求$p$为素数. 时间复杂度$O(\sqrt p \cdot \log(\sqrt{p}))$

int bsgs(int a,int b,int p){
	map<int,int>f;
	int m = ceil(sqrt(p));
	b %= p;
	for(int i = 1;i <= m;i++){
		b = 1ll*b*a%p;
		f[b] = i;
	}
	int v = Pow(a,m,p);
	b = 1;
	for(int i = 1;i <= m;i++){
		b = 1ll*b*v%p;
		if(f.count(b))return (1ll*i*m - f[b] + p)%p;
	}
	return -1;
}

求最小的正整数$n$使得$a^n \equiv 1 \pmod m$

int ord(int a,int m){
	a %= m;
	if(__gcd(a,m) != 1)return -1;
	int sq = sqrt(m)+1;
	map<int,int>mp;
	ll s = a;
	for(int i = 1;i <= sq;i++){
		if(!mp.count(s))mp[s] = i;
		s = s*a%m;
	}
	ll g = inv(Pow(a,sq,m),m),x = 1;
	for(int k = 0;k <= m/sq;k++){
		if(mp.count(x))return sq*k + mp[x];
		x = x*g%m;
	}
}

### 扩展版

$p$可以是任意模数. 前置:bsgs exgcd inv

int exbsgs(int a,int b,int p){
	if(b == 1 or p == 1)return 0;
	int g = __gcd(a,p),k = 0,na = 1;
	while(g > 1){
		if(b%g != 0)return -1;
		k++;
		b /= g;p /= g;na = na*(a/g)%p;
		if(na == b)return k;
		g = __gcd(a,p);
	}
	int f = bsgs(a,b*inv(na,p)%p,p);
	if(f == -1)return -1;
	return f+k;
}

## 拉格朗日插值

px,py为点坐标,p为模数,k是L(k)那个k inv为逆元

返回L(k)

ll lagrange(ll n,ll k){
	ll tans = 0;
	for(int i = 1;i <= n;i++){
		ll f = 1,g = 1;
		for(int j = 1;j <= n;j++)if(i != j){
			f = f*(k - px[j]+p)%p;
			g = g*(px[i] - px[j]+p)%p;
		}
		ll lota = py[i]*f%p*inv(g,p)%p;
		tans = (tans+lota)%p;
	}
	return tans; 
}

### 线性求$\sum_{i = 0}^{n}i^k$

(返回L(k)的值,多项式次数为n.按需更改yi的值,复杂度$O(k)$)

(注意求$\sum_{i = 0}^{n}i^k$需要调用line_lagrange(k+2,n) )

ll line_lagrange(ll n,ll k){
	ll tans = 0;
	get_inv();
	pre[0] = suf[n+1] = ifac[0] =  1;
	for(int i = 1;i <= n;i++){
		pre[i] = pre[i-1]*(k-i)%mod;
		ifac[i] = ifac[i-1]*invs[i]%mod;
	}
	for(int i = n;i >= 0;i--)suf[i] = suf[i+1]*(k-i)%mod;
	ll yi = 0; 
	for(int i = 1;i <= n;i++){
		yi = (yi+pow_mod(i,n-2))%mod;
		ll a = pre[i-1]*suf[i+1]%mod*yi%mod*ifac[i-1]%mod*ifac[n-i]%mod;
		if((n-i)%2)a = mod-a; 
		tans = (tans+a+mod)%mod;
	}
	return tans;
}

## 卢卡斯定理

Lucas(n,m) = C(n%p,m%p)*Lucas(n/p,m/p)%p

p需为素数,fac[]为预处理的阶乘,时间复杂度大约是$O(p + mlogm)$(预处理需要$O(p)$)

ll C(ll n,ll m){
	if(m > n)return 0;
	return (fac[n]*Pow(fac[m],p-2)%p*Pow(fac[n-m],p-2)%p)%p;
}
ll Lucas(ll n,ll m){
	if(!m)return 1;
	return C(n%p,m%p)*Lucas(n/p,m/p)%p;
}

线性基

maxbit取值域位数

insert返回0则说明现有基可以表示x

qmax返回异或最大值

struct LinerBasis{
	ll base[maxbit+4];
	bool vis[maxbit+4];
	bool insert(ll x){
		for(int i = maxbit;~i;i--){
			if(x&(1LL<<i)){
				if(!base[i]){
					base[i] = x;
					break;
				}
				else x ^= base[i];
			}
		}
		return x;
	}
	ll qmax(){
		ll ans = 0;
		for(int i = maxbit;~i;i--)ans = max(ans,ans^base[i]);
		return ans;
	}
}yay;

杜教BM(快速求模意义下线性递推式)

设$n$为输入项数, 复杂度$O(n^2\log n)$. 对于$k$阶递推式(例如斐波那契数列为2阶递推式)至少需要前$2k$项

#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,n) for (int i=a;i<n;i++)
#define pb push_back
typedef long long ll;
#define SZ(x) ((ll)(x).size())
typedef vector<ll> VI;
typedef pair<ll, ll> PII;
const ll mod = 1000000007;
ll powmod(ll a, ll b) { ll res = 1; a %= mod; assert(b >= 0); for (; b; b >>= 1) { if (b & 1)res = res * a % mod; a = a * a % mod; }return res; }
// head
 
ll _, n;
namespace linear_seq {
    const ll N = 10010;
    ll res[N], base[N], _c[N], _md[N];
 
    vector<ll> Md;
    void mul(ll* a, ll* b, ll k) {
        rep(i, 0, k + k) _c[i] = 0;
        rep(i, 0, k) if (a[i]) rep(j, 0, k) _c[i + j] = (_c[i + j] + a[i] * b[j]) % mod;
        for (ll i = k + k - 1; i >= k; i--) if (_c[i])
            rep(j, 0, SZ(Md)) _c[i - k + Md[j]] = (_c[i - k + Md[j]] - _c[i] * _md[Md[j]]) % mod;
        rep(i, 0, k) a[i] = _c[i];
    }
    ll solve(ll n, VI a, VI b) { // a 系数 b 初值 b[n+1]=a[0]*b[n]+...
//        printf("%d\n",SZ(b));
        ll ans = 0, pnt = 0;
        ll k = SZ(a);
        assert(SZ(a) == SZ(b));
        rep(i, 0, k) _md[k - 1 - i] = -a[i]; _md[k] = 1;
        Md.clear();
        rep(i, 0, k) if (_md[i] != 0) Md.push_back(i);
        rep(i, 0, k) res[i] = base[i] = 0;
        res[0] = 1;
        while ((1ll << pnt) <= n) pnt++;
        for (ll p = pnt; p >= 0; p--) {
            mul(res, res, k);
            if ((n >> p) & 1) {
                for (ll i = k - 1; i >= 0; i--) res[i + 1] = res[i]; res[0] = 0;
                rep(j, 0, SZ(Md)) res[Md[j]] = (res[Md[j]] - res[k] * _md[Md[j]]) % mod;
            }
        }
        rep(i, 0, k) ans = (ans + res[i] * b[i]) % mod;
        if (ans < 0) ans += mod;
        return ans;
    }
    VI BM(VI s) {
        VI C(1, 1), B(1, 1);
        ll L = 0, m = 1, b = 1;
        rep(n, 0, SZ(s)) {
            ll d = 0;
            rep(i, 0, L + 1) d = (d + (ll)C[i] * s[n - i]) % mod;
            if (d == 0) ++m;
            else if (2 * L <= n) {
                VI T = C;
                ll c = mod - d * powmod(b, mod - 2) % mod;
                while (SZ(C) < SZ(B) + m) C.pb(0);
                rep(i, 0, SZ(B)) C[i + m] = (C[i + m] + c * B[i]) % mod;
                L = n + 1 - L; B = T; b = d; m = 1;
            }
            else {
                ll c = mod - d * powmod(b, mod - 2) % mod;
                while (SZ(C) < SZ(B) + m) C.pb(0);
                rep(i, 0, SZ(B)) C[i + m] = (C[i + m] + c * B[i]) % mod;
                ++m;
            }
        }
        return C;
    }
    ll gao(VI a, ll n) {
        VI c = BM(a);
        c.erase(c.begin());
        rep(i, 0, SZ(c)) c[i] = (mod - c[i]) % mod;
        return solve(n, c, VI(a.begin(), a.begin() + SZ(c)));
    }
};
int main() {
    while (~scanf("%lld", &n)) {
        vector<ll>v;
        v.push_back(1);/*前几项*/v.push_back(5);v.push_back(15);
        v.push_back(35);v.push_back(70);v.push_back(126);
        v.push_back(210);v.push_back(330);v.push_back(495);v.push_back(715);
        v.push_back(1001);
        //输入n ,输出第n项的值  一般大于10项即可出答案，越多越好
        printf("%lld\n", linear_seq::gao(v, n - 1));
    }
}

特殊计数序列

第二类Stirling数

第二类Stirling数$S(n,m)$计数的是把$n$个元素的集合划分到$m$个不可区分的盒子中且没有空盒子的划分个数.

$S(n,n) = 1$, $S(n,0) = 0\space(n > 0)$

递推式: $S(n,m) = m \cdot S(n-1,m) + S(n-1,m-1)$

求第n行

\[ S(n,m) = \sum_{i = 0}^m \frac{i^n}{i!}\cdot\frac{(-1)^{m-i}}{(m-i)!} \] 用NTT求解即可. finv[i]为$i!$的逆元, 需要线性预处理出来. mul()见NTT板子. 时间复杂度$O(n \log n)$

Poly Stirling_row(int n){
	Poly A(n+1),B(n+1);
	for(int i = 0;i <= n;i++){
		A[i] = 1ll*Pow(i,n)*finv[i]%p;
		B[i] = ( ( 1ll * ((n-i)%2?-1:1) * finv[n-i] )%p+p)%p;
	}
	reverse(B.begin(),B.end());
	auto tans = mul(A,B);
	tans.resize(n+1);
	return tans;
}

求第k列

\[ \sum_{n = k}^{\infin}S(n,k) \cdot \frac{x^n}{n!} = \frac{(e^x - 1)^k}{k!} \]

因而有 \[ S(n,k) = \frac{n!}{k!}\cdot[n](e^x - 1)^k \] 多项式快速幂即可. 时间复杂度$O(n \log n)$

Poly Stirling_Col(int n,int k){
	Poly A(n+1);
	for(int i = 1;i <= n;i++)A[i] = finv[i];
	A = powP(A,k,k);
	A.resize(n+1);
	for(int i = 0;i <= n;i++)A[i] = 1ll*A[i]*finv[k]%p*fac[i]%p;
		return A;
}

Bell数

Bell数($B(n)$)计数的是将一个有$n$个元素的集合划分成任意个非空子集的方案数. 显然有$B(n) = \sum_{i = 1}^n S(n,i)$

前10个Bell数为1 2 5 15 52 203 877 4140 21147 115975 \[ B(n) = [n]e^{e^x - 1} \] 多项式指数函数即可. 时间复杂度$O(n\log n)$

Poly Bell(int n){
	Poly A(n+1);
	for(int i = 1;i <= n;i++)A[i] = finv[i];
	A = expP(A);
	for(int i = 0;i <= n;i++)A[i] = 1ll*A[i]*fac[i]%p;
	return A;
}

博弈论

威佐夫游戏

两人轮流取两堆筹码，其中取法有两种：取走一堆中任意个筹码，或从两堆中取走相同数目的筹码。取完所有筹码的一方获胜。

(1为先手胜)

if(n1>n2)  swap(n1,n2);
temp=floor((n2-n1)*(1+sqrt(5.0))/2.0);
if(temp==n1) cout<<"0"<<endl;
else cout<<"1"<<endl;

矩阵

定义&快速幂

注意：vector版的矩阵会比数组慢起码2倍，时间吃紧的时候慎用

数组版:

#define FOR(i,a,b) for(int i = a;i <= b;i++)
struct Matrix{
	ll A[maxn][maxn];
	void initer(){
        for(int i = 1;i <= n;i++)
            for(int j = 1;j <= n;j++)A[i][j] = 0;
		for(int i = 1;i <= n;i++)A[i][i] = 1;
	}
	Matrix operator * (const Matrix & a){
		Matrix res;
		memset(res.A,0,sizeof(res.A));
		FOR(k,1,n)
			FOR(i,1,n)
				FOR(j,1,n){
					res.A[i][j] = (res.A[i][j] + (A[i][k]*a.A[k][j]));
					if(res.A[i][j] > MOD2)res.A[i][j] -= MOD2;
				}
		FOR(i,1,n)
			FOR(j,1,n)res.A[i][j] %= MOD;
		return res;
	}
}

vector版:

typedef vector<vector<int> > Matrix;
void Minit(Matrix & a){
	int n = a.size();
	for(int i = 0;i < n;i++){
		for(int j = 0;j < n;j++)a[i][j] = (i == j);
	}
}
Matrix operator * (const Matrix & a,const Matrix & b){//需保证a的列数等于b的行数且两个矩阵都不为空
	int n = a.size(),m = b.size(),r = b[0].size();
	Matrix res(n,vector<int>(m));
	for(int k = 0;k < m;k++){
		for(int i = 0;i < n;i++){
			for(int j = 0;j < r;j++){
				res[i][j] = (res[i][j] + (1ll*a[i][k]*b[k][j])%MOD)%MOD;
			}
		}
	}
	return res;
}
Matrix matrixPow(Matrix x,ll d){
	int n = x.size();
	Matrix ans(n,vector<int>(n));
	Minit(ans);
	while(d){
		if(d&1)ans = ans*x;
		x = x*x;
		d /= 2;
	}
	return ans;
}

置换矩阵

每一行每一列有且仅有一个$1$, 其余元素为$0$的矩阵为置换矩阵.

置换矩阵的逆是它的转置.

循环矩阵

除第一行外, 其余每行都是上一行循环右移一格的矩阵是循环矩阵. 循环矩阵的乘积仍是循环矩阵.

因此只需保存和计算循环矩阵的第一行, 两个循环矩阵相乘复杂度为$O(n^2)$

typedef vector<int> cirMatrix;//循环矩阵只需保存第一行
Matrix buildCirMatrix(const cirMatrix & a){//构建循环矩阵,此处为右移
	int n = a.size();
	Matrix res(n,vector<int>(n));
	res[0] = a;
	for(int i = 1;i < n;i++){
		res[i][0] = res[i-1][n-1];
		for(int j = 1;j < n;j++)res[i][j] = res[i-1][j-1];
	}
	return res;
}
cirMatrix operator * (const cirMatrix & a,const cirMatrix & b){
	Matrix A = Matrix(1,a);
	Matrix B = buildCirMatrix(b);
	return (A*B)[0];
}
cirMatrix matrixPow(cirMatrix x,ll d){
	int n = x.size();
	cirMatrix ans(n);
	ans[0] = 1;
	while(d){
		if(d&1)ans = ans*x;
		x = x*x;
		d /= 2;
	}
	return ans;
}

高斯消元

输入一个$n$行$n+1$列的矩阵,其中$A[i][n]$代表等式右边第$i$个值. 复杂度$O(n^3)$

gauss(矩阵,n)返回1有解, 返回0无解

bool gauss(vector<vector<ld> > & A,int n){
	for(int i = 0;i < n;i++){
		int r = i;
		for(int j = i+1;j < n;j++)if(fabs(A[j][i]) > fabs(A[r][i]))r = j;
		if(r != i)for(int j = 0;j < n+1;j++)swap(A[i][j],A[r][j]);
		if(fabs(A[i][i]) <= eps)return 0;
		for(int j = n;j >= i;j--){
			for(int k = i+1;k < n;k++){//提高精度
				A[k][j] -= A[k][i]/A[i][i]*A[i][j];
			}
		}
	}
	for(int i = n-1;i >= 0;i--){
		for(int j = i+1;j < n;j++)A[i][n] -= A[j][n]*A[i][j];
		A[i][n] /= A[i][i];
	}
	return 1;
}

异或方程组

时间复杂度$O(\frac{n^3}{w})$, $w = 32或64$, 取决于评测机位数

const int maxb = 1005;
vector<bitset<maxb> > a(m,bitset<maxb>());
int xorRank(vector<bitset<maxb> > & a,int n,int m){
	//n个未知数, m个方程
	int mr = -1,i = 0,j = 0;
	while(i < m and j < n){
		int r = i;
		for(int k = i;k < m;k++)if(a[k][j]){r = k;break;}
		if(a[r].test(j)){
			mr = max(mr,r);
			if(r != i)swap(a[r],a[i]);
			for(int u = i+1;u < m;u++)if(a[u].test(j))a[u] ^= a[i];
			i++;
		}
		//else return 0; //某一行全为0, 说明有多解
		j++;
	}
	for(int k = m-1;k >= 0;k--){
		for(int j = k+1;j < m;j++)a[k][m-1] = (a[k][m-1] ^ (a[j][m-1] & a[k][j]));
	}
    return i;//返回该矩阵的秩
	//return mr+1;//返回最少使用了几个方程组来得到解
}

行列式

浮点数版

const ld eps = 1e-9;
ld caldet(vector<vector<ld> > & a,int n){
	ld det = 1;
	for(int i = 0;i < n;i++){
		int k = i;
		for(int j = i+1;j < n;j++){
			if(fabs(a[j][i]) > fabs(a[k][i]))k = j;
		}
		if(fabs(a[k][i]) < eps){
			det = 0;
			break;
		}
		swap(a[i],a[k]);
		if(i != k)det = -det;
		det *= a[i][i];
		for(int j = i+1;j < n;j++)a[i][j] /= a[i][i];
		for(int j = 0;j < n;j++){
			if(j != i and fabs(a[j][i]) > eps){
				for(int k = i+1;k < n;k++){
					a[j][k] -= a[i][k] * a[j][i];
				}
			}
		}
	}
	return det;
}

整数版

输入一个vector矩阵, 输出行列式模$p$的值. 均摊复杂度$O(n^3 + n^2\log n)$

int caldet(vector<vector<int> > & a,int n,int p){
	int det = 1,w = 1;
	for(int i = 0;i < n;i++){
		for(int j = i+1;j < n;j++){
			while(a[i][i]){
				int div = a[j][i]/a[i][i];
				for(int k = i;k < n;k++){
					a[j][k] = (a[j][k] - 1ll*div*a[i][k]%p + p)%p;
				}
				swap(a[i],a[j]);w = -w;
			}
			swap(a[i],a[j]);w = -w;
		}
	}
	for(int i = 0;i < n;i++)det = 1ll*a[i][i]*det % p;
	det = (1ll*w*det + p)%p;
	return det;
}

Kirchhoff 矩阵树定理

矩阵树定理用于求解图的生成树个数. 来源: https://oi-wiki.org/graph/matrix-tree/

设$t(G)$为图$G$的全部生成树边权积的和(当边权为1时即为全部生成树的个数)

即 \[ t(G) = \sum_{T \in G} \prod_{e \in T}w_e \] 其中$T$代表图$G$的生成树, $e$为生成树上的边，$w_e$为边的权值.

以下所讨论的图允许重边但不允许自环(去掉自环即可).

无向图形式-行列式

设$D$为无向图$G$的度数矩阵, $A$为无向图的邻接矩阵. $deg(x)$为点$x$的度数.

则$D[i][i] = deg(i) \cdot[i == j]$, $A[i][j] = |i与j间的连边权值之和|$

设$Laplace$矩阵$L = D - A$. $L_r$为$L$的$n-1$阶主子式(即将矩阵$L$去掉第$r$行第$r$)列.

则对任意的$r$, $L_r$的行列式为该无向图的生成树个数. 即$t(G) = det(L_r)$

无向图形式-特征值

若$\lambda_1,\lambda_2,...,\lambda_{n-1}$为$L$的$n-1$个非负特征值, 则有 \[ t(G) = \frac{1}{n}\prod_{i = 1}^{n-1}\lambda_i \]

有向图形式

有向图的生成树分两种: 根向, 所有边均由儿子指向父亲; 叶向, 所有边均由父亲指向儿子.

设$D^{out}$为有向图$G$的出度矩阵, $D^{in}$为有向图$G$的入度矩阵, 两者与无向图形式相同.

$A$为有向图邻接矩阵, $A[i][j] = |i指向j的边数|$

设入度Laplace矩阵$L^{in} = D^{in} - A$, 出度Laplace矩阵为$L^{out} = D^{out} - A$

则以$k$为根的叶向生成树个数为$t^{leaf}(G,k) = det(L_k^{in})$, 根向生成树个数为$t^{root}(G,k) =det(L_k^{out})$.

若要求所有生成树, 枚举$k$并求和即可.

带权形式

对有向图无向图均适用. 将带权边转化为无权的重边, 重边次数为边权即可.

BEST定理

通过图中所有边恰好一次且行遍所有顶点的回路称为欧拉回路.

具有欧拉回路的无向图或有向图称为欧拉图.

设 $G$ 是有向欧拉图，那么 $G$ 的不同欧拉回路总数 $ec(G)$ 是 \[ ec(G) = t^{root}(G,k)\prod_{v\in V}(\deg (v) - 1)! \]

注意，对欧拉图 $G$ 的任意两个节点 $k, k'$，都有 $t^{root}(G,k)=t^{root}(G,k')$，且欧拉图 $G$ 的所有节点的入度和出度相等。

一些例题

加速递推式

以$Fibonacci$为例:

$f(1) = f(2) = 1,f(n) = f(n-1) + f(n-2),(n \geq 10^9)$,求$f(n)$

设 \[ F_i = \begin{bmatrix} f_{i+1}\\ f_i \end{bmatrix}, F_0 = \begin{bmatrix} f_1\\ f_0 \end{bmatrix} \] 且$F_i = A^i \cdot F_0$

则有 \[ A = \begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix} \] 一般地, 对于递推式$f(i) = \sum_{i = 1}^{k}a_i \cdot f(i-k)$

有: \[ F_i = \begin{bmatrix} f(i+k-1)\\ f(i+k-2)\\ ...\\ f(i) \end{bmatrix}, A = \begin{bmatrix} a_1 & a_2 & a_3 & ... & a_k\\ 1 & 0 & 0& \cdots & 0\\ 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & 1 & 0 \end{bmatrix} \] 且$F_n = A^n \cdot F_0$

求包含严格k条边的方案数/最短路

1.无权图求从$i$到$j$经过$k$条边的方案数:求邻接矩阵的$k$次幂即可

2.有向图求包含$k$条边的最短路:求邻接矩阵的$k$次幂,把普通矩阵乘法改为Floyd

0-1分数规划

给出数列$\{a_i\},\{b_i\}$,求一组$w_i \in \{0,1\}$, 使得 \[ \frac{\sum_{i = 1}^na_i \cdot w_i}{\sum_{i = 1}^nb_i \cdot w_i} \] 最大或最小

或给出$n$个物品,每个物品有两个权值$a_i$和$b_i$, 选择任意个物品使得$\frac{\sum a_i}{\sum b_i}$最大

二分法

\[ \frac{\sum_{i = 1}^na_i \cdot w_i}{\sum_{i = 1}^nb_i \cdot w_i} \ge M \Leftrightarrow \sum_{i = 1}^n w_i \cdot (a_i - M\cdot b_i) \ge 0 \]

因此只需判断$_{i = 1}^n w_i (a_i - Mb_i) $的最大值是否非负即可对$M$二分.

Dinkelbach 算法(迭代法)

将上一轮的答案作为输入进行迭代.注意此时答案应为$\frac{\sum_{i = 1}^na_i}{\sum_{i = 1}^nb_i}$, 而不是变形后的式子.

例题

[POJ2728]Desert King

给出$n \leq 10^3$个点的横纵坐标和高度, 两点间距离为欧几里得距离, 价值为高度之差绝对值. 所有点构成一个完全图. 求一棵生成树使得

距离/价值比最小.

二分法

此题是有一定限制的0-1分数规划(所选的"物品"构成一棵生成树). 直接求边权为$w_i \cdot (a_i - M\cdot b_i)$的最小生成树, 判断是否有解即可二分.

由于是完全图, 采用不加优化的Prim求解.

复杂度大概是$O(n^2 \log((R - L)*2^{-eps}))$

由于二分法较慢, 本题需要猜一个答案上界(100)才能通过=-=

用时2610ms

//POJ2728
#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
typedef long long ll;
typedef double ld;
#define endl '\n'
const int maxn = 1e3+10;
const ld INF = 1e6+10;
const ld eps = 1e-5;
ld G[maxn][maxn],dis[maxn];
int x[maxn],y[maxn],h[maxn],n;
ld w[maxn][maxn];
bool vis[maxn];
ld Dis(int idx1,int idx2){
	return sqrt((x[idx1] - x[idx2])*(x[idx1] - x[idx2])+(y[idx1] - y[idx2])*(y[idx1] - y[idx2]));
}
bool prim(ld M){
	for(int i = 0;i <= n;i++)vis[i] = 0,dis[i] = INF;
	for(int i = 1;i < n;i++){
		int x = 0;
		for(int j = 1;j <= n;j++)if(!vis[j]){
			if(!x or dis[x] > dis[j])x = j;
		}
		vis[x] = 1;
		for(int j = 1;j <= n;j++)if(!vis[j]){
			dis[j] = min(dis[j],w[j][x] - M*G[j][x]);
		}
	}
	ld tans = 0;
	for(int i = 2;i <= n;i++)tans += dis[i];
	return tans >= 0.0;
}
signed main(){
	while(scanf("%d",&n) == 1){
		if(!n)break;
		for(int i = 1;i <= n;i++){
			scanf("%d%d%d",&x[i],&y[i],&h[i]);
		}
		for(int i = 1;i <= n;i++){
			for(int j = 1;j <= n;j++){
				G[i][j] = Dis(i,j);
				w[i][j] = abs(h[i] - h[j]);
			}
		}
		ld M = 0,L = 0,R = 100;
		while(R - L >= eps){
			M = (L+R)/2;
			if(prim(M))L = M;
			else R = M;
		}
		printf("%.3f\n",L);
	}
	return 0;
}

迭代法

求生成树时仍按边权$w_i \cdot (a_i - M\cdot b_i)$来计算. 但统计答案则计算$\frac{\sum_{i = 1}^na_i}{\sum_{i = 1}^nb_i}$

跑的飞快(250ms)

//POJ2728
#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
typedef long long ll;
typedef double ld;
#define endl '\n'
const int maxn = 1e3+10;
const ld INF = 1e9+10;
const ld eps = 1e-5;
ld G[maxn][maxn],dis[maxn];
int x[maxn],y[maxn],h[maxn],to[maxn],n;
ld w[maxn][maxn];
bool vis[maxn];
ld Dis(int idx1,int idx2){
	return sqrt((x[idx1] - x[idx2])*(x[idx1] - x[idx2])+(y[idx1] - y[idx2])*(y[idx1] - y[idx2]));
}
ld prim(ld M){
	ld up = 0,down = 0;
	for(int i = 0;i <= n;i++)vis[i] = 0,dis[i] = INF,to[i] = 1;
	for(int i = 1;i < n;i++){
		int x = 0;
		for(int j = 1;j <= n;j++)if(!vis[j]){
			if(!x or dis[x] > dis[j])x = j;
		}
		vis[x] = 1;
		for(int j = 1;j <= n;j++)if(!vis[j]){
			if(dis[j] > w[x][j] - M*G[x][j])dis[j] = w[x][j] - M*G[x][j],to[j] = x;
		}
	}
	for(int i = 2;i <= n;i++){
		up += w[i][to[i]];
		down += G[i][to[i]];
	}
	return up/down;
}
signed main(){
	while(scanf("%d",&n) == 1){
		if(!n)break;
		for(int i = 1;i <= n;i++){
			scanf("%d%d%d",&x[i],&y[i],&h[i]);
		}
		for(int i = 1;i <= n;i++){
			for(int j = 1;j <= n;j++){
				G[i][j] = Dis(i,j);
				w[i][j] = abs(h[i] - h[j]);
			}
		}
		ld ans = 0,lans = 1;
		while(fabs(lans - ans) > eps){
			lans = ans;
			ans = prim(ans);
		}
		printf("%.3f\n",ans);
	}
	return 0;
}

多项式

FFT

复数类操作:

1
2
3

complex<double>m;
m.real();m.imag();//获取实部虚部的值
m.real(114);m.imag(514);//为实部虚部赋值

包括手写复数类, DFT与IDFT, 多项式乘法.

复数类: $x$为实部, $y$为虚部
FFT(Comp * A,int siz,int type)

A: 复数数组, 全局变量, 为待求解的多项式

siz: 多项式长度. 需保证为2的次方

type: 1为DFT, -1为IDFT
Poly mul(Poly A,Poly B)

多项式乘法. 输入两个多项式$A,B$, 长度分别为$n,m$(也就是$n-1,m-1$阶). 返回一个长度$n - m + 1$的多项式, 为$A,B$的乘积.

Poly为vector<int>.

!!!切记在点值表示法下相乘时, 两个IDFT的长度必须相同!!!

typedef double ld;
typedef vector<int> Poly;
const ld PI = acos(-1.0);
int rev[maxn];
struct Comp{
    double x, y;
    Comp (double x = 0, double y = 0):x(x),y(y){}
    Comp operator * (Comp B){
   		return Comp(x*B.x - y*B.y,x*B.y + y*B.x);
	}
	Comp operator + (Comp B){
	    return Comp(x + B.x,y + B.y);
	}
	Comp operator - (Comp B){
	    return Comp(x - B.x,y - B.y);
	}
};
Comp F1[maxn],F2[maxn];
void FFT(Comp * A,int siz,int type){
	int n = siz;
	int S = log2(n);
	for(int i = 0;i < n;i++)rev[i] = (rev[i>>1] >> 1) | ((i & 1) << (S - 1));
	for(int i = 0;i < n;i++)if(i < rev[i])swap(A[i],A[rev[i]]);
	for(int i = 1;i < n;i *= 2){
		Comp Wn(cos(PI/i),type*sin(PI/i));
		for(int j = 0;j < n;j += i*2){
			Comp W(1,0);
			for(int k = 0;k < i;k++){
				Comp facx = A[j+k],facy = W*A[j+k+i];
				A[j+k] = facx + facy;
				A[j+k+i] = facx - facy;
				W = W*Wn;
			}
		}
	}
    if(type == -1)for(int i = 0;i < n;i++)A[i].x = (int)((A[i].x/n + 0.5));
}
Poly mul(Poly A,Poly B){
	int n = A.size(),m = B.size();
	int siz = n + m - 1;
	Poly C(siz);
	if(siz < 64){//小于等于64项直接暴力算
		for(int i = 0;i < n;i++){
			for(int j = 0;j < m;j++)C[i+j] = (C[i+j] + 1LL*A[i]*B[j]);
		}
		return C;
	}
	int fsiz = 1;
	while(fsiz <= siz)fsiz *= 2;
	for(int i = 0;i < fsiz;i++)F1[i] = F2[i] = 0;
	for(int i = 0;i < n;i++)F1[i] = A[i];
	for(int i = 0;i < m;i++)F2[i] = B[i];
	FFT(F1,fsiz,1);
	FFT(F2,fsiz,1);
	for(int i = 0;i < fsiz;i++)F1[i] = F1[i]*F2[i];
	FFT(F1,fsiz,-1);
	for(int i = 0;i < siz;i++){
		C[i] = ((ll)F1[i].x);
	}
	return C;
}

卷积

\[ C[i] = \sum_{j = 0}^iA[j]B[i-j] \]

卷积等同于多项式乘积

NTT

从OI-wiki抄来的, 用法和上面的FFT一致. 记得自己手写一遍

const int maxn = 6e6+10;
const int p = 998244353;
int Pow(int x,int d){
	int tans = 1;
	if(d == 0)return 1%p;
	int a = Pow(x,d/2);
	tans = 1ll*a*a%p;
	if(d%2)tans = 1ll*tans*x%p; 
	return tans%p;
}
typedef vector<int> Poly;
int F1[maxn],F2[maxn];
int rev[maxn];
void NTT(int * A,int lim,int opt) {
 	int i, j, k, m, gn, g, tmp;
 	for(int i = 0; i < lim; ++i)rev[i] = (i & 1)*(lim >> 1) + (rev[i >> 1] >> 1);
 	for(i = 0; i < lim; ++i)if (rev[i] < i) swap(A[i], A[rev[i]]);
 	for(m = 2; m <= lim; m <<= 1) {
	    k = m >> 1;
	    gn = Pow(3,(p - 1) / m);
	    for(i = 0; i < lim; i += m) {
	    	g = 1;
	    	for (j = 0; j < k; ++j, g = 1ll * g * gn % p) {
	    	    tmp = 1ll * A[i + j + k] * g % p;
	    	    A[i + j + k] = (A[i + j] - tmp + p) % p;
	     	 	A[i + j] = (A[i + j] + tmp) % p;
	    	}
		}
	}
	if(opt == -1){
    	reverse(A+1,A+lim);
    	int inv = Pow(lim,p-2);
    	for(i = 0; i < lim; ++i) A[i] = 1ll * A[i] * inv % p;
	}
}
Poly mul(Poly A,Poly B){
	int n = A.size(),m = B.size();
	int siz = n + m - 1;
	Poly C(siz);
	if(siz < 64){//小于等于64项直接暴力算
		for(int i = 0;i < n;i++){
			for(int j = 0;j < m;j++)C[i+j] = (C[i+j] + 1LL*A[i]*B[j]%p)%p;
		}
		return C;
	}
	int fsiz = 1;
	while(fsiz <= siz)fsiz *= 2;
	for(int i = 0;i < fsiz;i++)F1[i] = F2[i] = 0;
	for(int i = 0;i < n;i++)F1[i] = A[i];
	for(int i = 0;i < m;i++)F2[i] = B[i];
	NTT(F1,fsiz,1);
	NTT(F2,fsiz,1);
	for(int i = 0;i < fsiz;i++)F1[i] = 1ll*F1[i]*F2[i]%p;
	NTT(F1,fsiz,-1);
	for(int i = 0;i < siz;i++){
		C[i] = F1[i];
	}
	return C;
}

与FFT的区别在点值表示法相乘时, NTT需要取模. 另外注意原根需与模数对应

1	for(int i = 0;i < (int)A.size();i++)A[i] = 1llA[i]B[i]%p;

分治FFT/NTT

设有$n$个多项式, 第$i$个为$Func[i]$

就嗯分治完事了. 时间复杂度$O(n\log^2 n)$.

vector<Poly>Func;
Poly cal(int L,int R){
	if(R - 1 == L)return Func[L];
	int M = (L+R)/2;
	return mul(cal(L,M),cal(M,R));
}

求逆

板子来源: Siyuan -「算法笔记」多项式求逆

调用invP(Poly A),返回模$x^n$意义下$A$的逆. 其中$A$是一个$n-1$次多项式(实际上是一个长度为$n$的vector<int>)

时间复杂度: $O(nlogn)$

int P[maxn],Q[maxn];
Poly invP(Poly A){
	int n = A.size();
	int siz = 1;
	while(siz < n)siz *= 2;
	A.resize(siz);
	Poly I(siz);
	I[0] = Pow(A[0],p-2);//inv(A[0]);
	for(int h = 2;h <= siz;h *= 2){
		int H = h*2;
		for(int i = 0;i <= H;i++)P[i] = Q[i] = 0;//重要,必须清空
		for(int i = 0;i < h;i++)P[i] = A[i],Q[i] = I[i];
		NTT(P,H,1);
		NTT(Q,H,1);
		for(int i = 0;i < H;i++){
			P[i] = 1ll*Q[i]*(2 - 1ll*P[i]*Q[i]%p + p)%p;
		}
		NTT(P,H,-1);
		for(int i = 0;i < h;i++)I[i] = P[i];
	}
	I.resize(n);
	return I;
}

求导, 积分

derP求导, intP积分, 复杂度均为$O(n)$. 其中积分需要先预处理逆元inv[]

Poly derP(Poly A,bool samelen = 1){
	int n = A.size();
	Poly D(n-1);
	for(int i = 1;i < n;i++)D[i-1] = 1ll*i*A[i]%p;
	if(samelen)D.resize(n);
	return D;
}
Poly intP(Poly A,bool samelen = 1){
	int n = A.size();
	Poly ans(n+1);
	for(int i = 1;i <= n;i++)ans[i] = 1ll*inv[i]*A[i-1]%p;
	if(samelen)ans.resize(n);
	return ans;
}

取对数

输入$A(x)$, 求$B(x) \equiv ln(A(x)) \pmod {x^n}$. 两边同时求导得: \[ B'(x) = \frac{A'(x)}{A(x)}\text dx \] 再积分: \[ B(x) = ln(A(x)) = \int B'(x)\text dx = \int\frac{A'(x)}{A(x)}\text dx \] 需要求逆,求导和积分, 需保证$A[0] = 1$. 复杂度$O(n\log n)$

Poly lnP(Poly A){
	Poly tans = intP(mul(derP(A),invP(A)));
	tans.resize((int)A.size());
	return tans;
}

指数函数

前置: 求$\ln$

和求逆基本相同, 返回$e^{A(x)}$. 复杂度$O(nlogn)$

Poly expP(Poly A){
	int n = A.size();
	int siz = 1;
	while(siz < n)siz *= 2;
	A.resize(siz);
	Poly E(siz);
	E[0] = 1;
	for(int h = 2;h <= siz;h *= 2){
		Poly B(h),C(h);;
		for(int i = 0;i < h;i++)B[i] = C[i] = E[i];
		B = lnP(B);
		for(int i = 0;i < h;i++)B[i] = (1ll*-B[i] + A[i] + 2ll*p)%p;
		B[0] = (B[0] + 1)%p;
		B = mul(B,C);
		for(int i = 0;i < h;i++)E[i] = B[i];
	}
	E.resize(n);
	return E;
}

快速幂

求$A(x)^k$, 显然有$A(x)^k = e^{k\cdot ln(A(x))}$

保证$A[0] = 1$时:

Poly powP(Poly A,int k){
	A = lnP(A);
	for(int i = 0;i < (int)A.size();i++)A[i] = 1ll*A[i]*k%p;
	return expP(A);
}

$A[0] \neq 1$时:

首先判断是否有$k \ge n$且$A[0] == 0$. 若是, 则返回一个全为$0$的多项式. 注意, 必须在取模前判断大小, 且下面的模板没有包含这部分内容

然后对指数取模, 若输入的指数$k > p$, 则令k1 = k % p, k2 = k % (p - 1) (即$k2$对$\varphi(p)$取模)

Poly powP(Poly A,int k1,int k2){
	int n = A.size(),d = 0;
	while(d < n and !A[d])d++;
	if(1ll*d*k1 >= n)return Poly(n);
	A.erase(A.begin(),A.begin()+d);
	int k = Pow(A[0],p-2);
	for(int i = 0;i < (int)A.size();i++)A[i] = 1ll*A[i]*k%p;
	A = powP(A,k1);
	A.resize(n);
	k = Pow(k,p - 1 - k2);
	for(int i = 0;i < (int)A.size();i++)A[i] = 1ll*A[i]*k%p;
	d *= k1;
	for(int i = n-1;i >= 0;i--){
		A[i] = (i >= d?A[i-d]:0);
	}	
	return A;
}

高阶差分/前缀和

求$A[]$的$k$阶差分/前缀和. 设$A$的OGF: $F(x) = \sum_{i = 0}^{\infin}A_ix^i$

显然前缀和是$F(x)\cdot\sum_{i=0}^{\infin}x^i = F(x)\cdot \frac{1}{1-x}$, 差分是$F(x)\cdot(1-x)$

那么$k$阶前缀和是$F(x)\cdot (\frac{1}{1-x})^k$, 差分是$F(x)\cdot(1-x)^k$. 复杂度为$O(nlogk)$, 但常数巨大且很不好写.

另有结论:

前缀和$F(x)\cdot\sum_{i = 0}^nC_{k+i-1}^ix^i$

差分$F(x)\cdot\sum_{i=0}^n(-1)^iC_n^ix^i$

复杂度$O(n\log n)$. 使用前先预处理逆元.

Poly presumP(Poly A,ll k){
	int n = A.size();
	Poly G(n);
	G[0] = 1;
	for(int i = 1;i < n;i++)G[i] = G[i-1]*(k+i-1)%p*inv[i]%p;
	A = mul(A,G);
	return A;
}
Poly diffP(Poly A,ll k){
	int n = A.size();
	Poly G(n);
	G[0] = 1;
	for(int i = 1;i < n;i++)G[i] = (G[i-1]*(-1) + p)%p*(k-i+1)%p*inv[i]%p;	
	A = mul(A,G);
	return A;
}

前缀和与差分

区间加多项式

给出$f(x)$和$[L,R]$, 对$a[L...R]$分别加上$f(1),f(2),...,f(R-L+1)$.

结论: $k$阶多项式的$k+1$阶差分余项为常数. 对$k+1$阶差分后的余项做$k+1$次前缀和即可还原原多项式.

为了维护该操作, 我们先对原数组做$k+1$阶差分, 然后算出$f(1),...,f(k+1)$的$k+1$阶差分. 此后可以按一般的差分来维护.

细节见代码

int d[maxn];
void diff(vector<int>&a,int cnt){
	while(cnt--)for(int i = (int)a.size() - 1;i > 0;i--){a[i] -=a[i-1];if(a[i] < 0)a[i] += p;}
}
void pre(vector<int>&a,int cnt){
	while(cnt--)for(int i = 1;i < (int)a.size();i++)a[i] = (a[i] + a[i-1])%p;
}
ll f(ll x,vector<ll>fac){
	ll res = 0,times = 1;
	for(ll v:fac){
		res = (res + v*times%p)%p;
		times = times*x%p;
	}
	return res;
}
ll g(ll x,vector<ll>fac,ll len){
	return (-f(x+len,fac) + p)%p;
}
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.setf(ios::fixed,ios::floatfield);
	cout.precision(12);	 
	int n,m,q;
	cin >> n >> m >> q;
	vector<int>a(n+100);
	int num = 6;
	for(int i = 0;i < n;i++)cin >> a[i];
	diff(a,num);
	while(m--){
		int l,r,k;
		cin >> l >> r >> k;
		l--,r--;
		vector<ll>fac(k+1);
		for(int i = k;i >= 0;i--)cin >> fac[i];
		vector<int>add(10+1),sub(10+1);
		for(int i = 0;i <= 10;i++)add[i] = f(i+1,fac),sub[i] = g(i+1,fac,r - l + 1);
		diff(add,num);
		diff(sub,num);
		for(int i = 0;i <= 10;i++)a[i+l] += add[i],a[i+r+1] += sub[i],a[i+l] %= p,a[i+r+1] %= p;
	}
	pre(a,num+1);
	while(q--){
		int l,r;
		cin >> l >> r;
		l--,r--;
		ll ans = (a[r] - (l?a[l-1]:0) + p)%p;
		cout << ans << endl;
	}
	return 0;
}

组合数学

杂乱的结论

长度为n, 可选[1,m]的本质不同序列个数

(本质不同即将序列排序后两个序列不相同) \[ C_{n+m-1}^{m-1} \]

递推组合数

\[ C_n^0 = 1\\ C_n^i = C_n^{i-1}\cdot \frac{n-i+1}{i}\\ C_{n+i-1}^i = C_{n+i-2}^{i-1}\cdot \frac{n+i-1}{i} \]

常用公式

$$ \[\begin{align} \sum_{i = 0}^n C_n^i &= 2^n\\ \sum_{i = 0}^n C_n^i \cdot x^i &= (x+1)^n(二项式定理)\\ \sum_{i = 0,i += 2}^n C_n^i &= 2^{n-1}\\ \sum_{i=1}^nC_n^i \cdot i &= n \cdot 2^{n-1}\\ \sum_{i = 0}^r C_m^i \cdot C_n^{r - i} &= C_{n+m}^r\\ \sum_{i = 0}^n (C_n^i)^2 = C_{2n}^n\\ \end{align}\] $$

常用结论

\[ \sum_{i = 1}^ni*m^{n-i} *C(n,i) = n(m+1)^{n-1} \]

生成函数

OGF

\[ f(x) = \sum_{i = 0}^{\infin}h_ix^i \]

常用公式

Part1

\[ \begin{align} \sum_{i=0}^nx^i &= \frac{1-x^{n+1}}{1-x} (等比数列求和)\\ \sum_{i=0}^\infin x^i &= \frac{1}{1-x}\\ \sum_{i = 0}^\infin x^{c\cdot i} &= \frac{1}{1 - x^c}\\ \sum_{i = 0}^{\infin}(i+1)\cdot x^i &= (\frac{1}{1-x})^2 \end{align} \]

Part2 \[ \frac{1}{(1-x)^n} = \sum_{i = 0}^\infin C_{n+i-1}^{i-1}x^i\\ \frac{1}{(1-x)^{n+1}} = \sum_{i = 0}^\infin C_{n+i}^i x^i \]

EGF

\[ g(x) = \sum_{i=0}^{\infin}h_x\cdot\frac{x^i}{i!} \]

常用公式

\[ \begin{align} \sum_{i = 0}^{\infin}\frac{x^i}{i!} & = e^x\\ \sum_{i = 0}^{\infin} (-1)^i \frac{x^i}{i!} & = e^{-x} \\ \sum_{i = 0,i +=2}^{\infin}\frac{x^i}{i!} &= \frac{e^x + e^{-x}}{2}\\ \sum_{i = 1,i += 2}^{\infin}\frac{x^i}{i!} & = \frac{e^x - e^{-z}}{2}\\ \sum_{i = 1}^{\infin}(-1)^{i+1} \frac{x^i}{i} & = \ln{(1+x)} \end{align} \]

反演

整除分块

设$f(x) = \lfloor \frac{k}{x} \rfloor$

$f(x)$分布呈块状,对于任意一个$i$,有最大的$j = \lfloor \frac{k}{\lfloor \frac{k}{i} \rfloor} \rfloor$,使得$f(i) == f(i+1) == ... = f(j)$

对于类似$\sum_{i = 1}^{n}\lfloor \frac{k}{i} \rfloor$的式子,可分块在$O(\sqrt{n})$时间内计算完毕

代码:

for(int l = 1;l <= n;l = r+1){
		if(k/l)r = min(n,k/(k/l));
		else r = n;
		ans += (r-l+1)*(k/l);
	}

对于向上取整, 有 \[ \lceil \frac{k}{x} \rceil = \lfloor \frac{k - 1}{x}\rfloor + 1 \]

狄利克雷(Dirichlet)卷积

定义

设$f,g$为两个数论函数,其狄利克雷卷积$f*g$为: \[ f \ast g(n) = \sum_{d|n}f(d)g(\frac{n}{d}) \]

常用式子

$$ \[\begin{aligned} \varepsilon=\mu \ast 1&\iff\varepsilon(n)=\sum_{d\mid n}\mu(d)\\ d=1 \ast 1&\iff d(n)=\sum_{d\mid n}1\\ \sigma=\text{id} \ast 1&\iff\sigma(n)=\sum_{d\mid n}d\\ \varphi=\mu \ast \text{id}&\iff\varphi(n)=\sum_{d\mid n}d\cdot\mu(\frac{n}{d}) \\ \end{aligned}\]

常用结论

\[ \sum_{d|n}^n\varphi(d) = n \]

\[ [\gcd(i,j) == 1] = \epsilon(\gcd(i,j) = \sum_{d|\gcd(i,j)}\mu(d) \]

套路题

$\sum_{i=1}^ngcd(i,n)$

\[ \begin{align} &\sum_{i=1}^n\gcd(i,n)\\ &=\sum_{d|n}d\sum_{i=1}^n[\gcd(i,n)==d]\\ &=\sum_{d|n}d\sum_{i=1}^{\frac{n}{d}}[\gcd(i,n)==1]\\ &=\sum_{d|n}d\cdotφ(\frac{n}{d}) \end{align} \]

$\sum_{i = 1}^n\sum_{j=1}^n\gcd(i,j)$

\[ \begin{align} &\sum_{i = 1}^n\sum_{j=1}^n\gcd(i,j)\\ &= \sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j) == d]\\ &=\sum_{d=1}^nd\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[gcd(i,j)==1]\\ &=\sum_{d=1}^nd(2\cdot\sum_{i=1}^{\frac{n}{d}}\varphi(i) - 1) \end{align} \]

关于最后一步,我们设$f(x) = \sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)==1]$

手画一下,容易看出$f(x) = 2\cdot\sum_{i=1}^n\varphi(i) - 1$

另一种做法:

我们知道$\varphi$有这么一个性质,$\sum_{d|n}^{n}\varphi(d) = n$

也就是$\sum_{d|\gcd(i,j)}^{n}\varphi(d) = \gcd(i,j)$

那么: \[ \begin{align} &\sum_{i = 1}^n\sum_{j = 1}^n\gcd(i,j)\\ &=\sum_{i=1}^n\sum_{j=1}^n\sum_{d|\gcd(i,j)}\varphi(d)\\ &=\sum_{d = 1}^n\varphi(d)\sum_{i=1}^n[d|i]\sum_{j=1}^n[d|j]\\ &=\sum_{d = 1}^n\varphi(d)\cdot\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor \end{align} \] 少处理个前缀和=-=而且我感觉这个做法更好理解

例题:洛谷P2398

$\sum_{i = 1}^{n}\sum_{j=i+1}^{n}\gcd(i,j)$

和上面那个基本一样 \[ \begin{align} &\sum_{i = 1}^n\sum_{j=i+1}^n\gcd(i,j)\\ &= \sum_{d=1}^nd\sum_{i=1}^n\sum_{j=i+1}^n[gcd(i,j) == d]\\ &=\sum_{d=1}^nd\sum_{i=1}^{\frac{n}{d}}\sum_{j=i+1}^{\frac{n}{d}}[gcd(i,j)==1]\\ &=\sum_{d=1}^nd(\sum_{i=1}^{\frac{n}{d}}\varphi(i) - 1) \end{align} \] (最后一步的化简和上一题大同小异,手画一下就有了)

$\sum_{i=1}^{n}\sum_{i=1}^{m}[\gcd(i,j) == 1]$

前置知识: \[ \epsilon(n)=\sum_{d|n}\mu(d)= \begin{cases} 1&n=1\\ 0&otherwise \end{cases} \]

开始愉快地推式子: \[ \begin{align} &\sum_{i=1}^{n}\sum_{i=1}^{m}[\gcd(i,j) == 1] \\ &=\sum_{i=1}^{n}\sum_{i=1}^{m}\epsilon(\gcd(i,j))\\ &=\sum_{i=1}^{n}\sum_{i=1}^{m}\sum_{d|gcd(i,j)}\mu(d)\\ &=\sum_{d = 1}^{\min(n,m)}\mu(d)\sum_{i = 1}^n[d|i]\sum_{j = 1}^{m}[d|j]\\ &=\sum_{d=1}^{\min(n,m)}\mu(d)\lfloor \frac{n}{d}\rfloor \lfloor \frac{m}{d} \rfloor \end{align} \] 同理,我们有: \[ \sum_{i = 1}^n\sum_{j = 1}^n\sum_{k = 1}^n [\gcd(i,j,k) == 1] = \sum_{d = 1}^n \mu(d)(\lfloor \frac{n}{d}\rfloor)^3 \]

杂项

小学数学

\[ \sum a_i \cdot b_j = \frac{\sum a_i^2 - \sum b_j^2}{2} \]

因为$(a+b)^2 = a^2 + b^2 + 2ab$

\[ \sum_{i = 1}^n i = \frac{n \cdot(n+1)}{2}\\ \sum_{i = 1}^n i^2 = \frac{n \cdot (n+1) \cdot(2n + 1)}{6}\\ \sum_{i = 1}^n i^3 = (\frac{n \cdot (n+1)}{2})^2 \]

杂乱的知识点

\[ \lfloor \frac{a}{bc} \rfloor = \lfloor \frac{\lfloor \frac{a}{b} \rfloor}{c} \rfloor \]
\[ a +b = a\oplus b + 2\cdot a \&b = a|b + a \&b \]
设$c(n)$为第$n$个Catalan数, 则有$\sum_{k = 0}^{n-1}c(k)\cdot c(n-1-k) = c(n)$
区间加平方和(范围$[l,n]$)可通过差分进行. 设差分数组为$d$, 每次对$[l,n]$加平方和的操作转换为d[l]++,d[l+1]++, 然后求d的3阶前缀和即可单点查询.
对于一个$n$个点的竞赛图(有向完全图), 不经过重复边的最长路最大为$(n为偶数): \frac{n\cdot (n-1)}{2} - \frac{n}{2} + 1$;$(n为奇数):\frac{n\cdot (n-1)}{2}$; 最小为$n-1$

图兰定理

对于一个有 $n$个点的无向图，若其中不存在三个点的环，则边数不超过 $n^2/4$

Pick定理

给定一个顶点均为整点(坐标为整数的点)的简单多边形,其面积$A$和内部格点数$I$,边上格点数$B$的关系是: $A = I+B/2 - 1$

Cayley公式

一个有$n$个节点的完全图有$n^{n-2}$种不同的生成树

Polya定理1

-仅旋转,染色数和顶点数均为n.欧拉函数优化(抄自Lskkkno1的题解qwq)

ll polya(int n,int k){
	ll tans = 0;
	for(int i = 1;i*i <= n;i++){
		if(n%i)continue;
		tans += (cal_phi(n/i)*pow_mod(n,i-1))%MOD;
		if(i*i != n)tans += cal_phi(i)*pow_mod(n,n/i-1)%MOD;
	}
	return tans%MOD;
}

\[ C_k^k +C_{k+1}^k +...+C_{n}^k（n > k) \]

某些数列

对于数列$\{1,2,2,3,3,3,4,4,4,4,...\}$(数$i$出现了$i$)次, 有通项公式 \[ a_n = \frac{1+2\sqrt{2n}}{2} \]

\[ \sum_{i = 1}^\infin i \cdot \frac{1}{a^i} \]

海伦公式

$a,b,c$为三角形边长, $s = \frac{a+b+c}{2}$

三角形面积为$S = \sqrt{s(s-a)(s-b)(s-c)}$

组合数奇偶性

对于$C_n^k$, 如果$n\&k == k$为奇数, 否则偶数.(按位and)

生成2进制恰好包含k位的下一个数

若$x$里有$k$个1, 则gen_next(x)为包含$k$个1的下一位数

int gen_next(int x){
	int b = x & -x;
	int t = x + b;
	int c = t & -t;
	int m = (c/b >> 1) - 1;
	return t | m;
}

求区间[n,m]内不含平方因子的数

筛掉$p^2$即可. $p$是素数且$p \in [n,m]$

等比数列

$a$为首项, $q$为公比

\[ 求和:S_n = a \cdot \frac{1-q^n}{1-q} (q \ne 1) \]

\[ 求积:P_n = a^n \cdot q^{\frac{n(n-1)}{2}} \]

大素数分解

namespace prime_fac {

    const int S = 8; // 随机算法判定次数，8~10 就够了

// 龟速乘
    long long mult_mod(long long a, long long b, long long c) {
        a %= c, b %= c;
        long long ret = 0;
        long long tmp = a;
        while (b) {
            if (b & 1) {
                ret += tmp;
                if (ret > c) ret -= c;
            }
            tmp <<= 1;
            if (tmp > c) tmp -= c;
            b >>= 1;
        }
        return ret;
    }

    // 快速幂
    long long qow_mod(long long a, long long n, long long _mod) {
        long long ret = 1;
        long long temp = a % _mod;
        while (n) {
            if (n & 1) ret = mult_mod(ret, temp, _mod);
            temp = mult_mod(temp, temp, _mod);
            n >>= 1;
        }
        return ret;
    }

    // 是合数返回true，不一定是合数返回false
    bool check(long long a, long long n, long long x, long long t) {
        long long ret = qow_mod(a, x, n);
        long long last = ret;
        for (int i = 1; i <= t; i++) {
            ret = mult_mod(ret, ret, n);
            if (ret == 1 && last != 1 && last != n - 1) return true;
            last = ret;
        }
        if (ret != 1) return true;
        return false;
    }

    // 是素数返回true，不是返回false
    mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
    bool Miller_Rabin(long long n) {
        if (n < 2) return false;
        if (n == 2) return true;
        if ((n & 1) == 0) return false;
        long long x = n - 1;
        long long t = 0;
        while ((x & 1) == 0) { x >>= 1; t++; }

        for (int i = 0; i < S; i++) {
            long long a = rng() % (n - 1) + 1;
            if (check(a, n, x, t))
                return false;
        }

        return true;
    }

    long long factor[100];// 存质因数
    int tol; // 质因数的个数，0~tol-1

    long long gcd(long long a, long long b) {
        long long t;
        while (b) {
            t = a;
            a = b;
            b = t % b;
        }
        if (a >= 0) return a;
        return -a;
    }

    long long pollard_rho(long long x, long long c) {
        long long i = 1, k = 2;
        long long x0 = rng() % (x - 1) + 1;
        long long y = x0;
        while (1) {
            i++;
            x0 = (mult_mod(x0, x0, x) + c) % x;
            long long d = gcd(y - x0, x);
            if (d != 1 && d != x) return d;
            if (y == x0) return x;
            if (i == k) { y = x0; k += k; }
        }
    }
    // 对n质因数分解，存入factor，k一般设置为107左右
    void findfac(long long n, int k) {
        if (n == 1) return;
        if (Miller_Rabin(n)) {
            factor[tol++] = n;
            return;
        }
        long long p = n;
        int c = k;
        while (p >= n) p = pollard_rho(p, c--);
        findfac(p, k);
        findfac(n / p, k);
    }
    vector<int> fac(long long n) {
        tol = 0;
        vector<int>ret;
        findfac(n, 107);
        for (int i = 0; i < tol; i++)ret.push_back(factor[i]);
        return ret;
    }
}

子集

子集的和与积

对于集合$S = \{a_i\}，i = 1...n$, 其全部子集中元素和的和为$2^{n-1}\cdot \sum_{i = 1}^na_i$
其全部子集中元素积的和为$\prod_{i = 1}^n(a_i + 1)$
从中选出大小为$k$的全部子集, 元素积的和可用OGF求解, 即 \[ \sum_{S' \in S \and |S'| = k}\space \prod_{p_i\in S'}p_i = [k] \prod_{i = 1}^n (1+a_i) \] NTT分治即可.

枚举子集

给定集合$m$, 枚举它的全部子集:
1
for(int S = m;S >= 0;S = (S-1)&m)
枚举大小为$n$的全部集合$m$, 且对于每个$m$, 枚举其全部子集$S$. 时间复杂度$O(3^n)$:
1
2
3
4
for(int m = 0;m < (1<<n);m++){
for(int S = m;S >= 0; S = (S-1)&m)
}

子集和/高维前缀和

求解如下问题:

子集和: \[ f[S] = \sum_{T \in S}f[T] \] 超集和: \[ f[T] = \sum_{T\in S}f[S] \]

for(int i = 0; i < n; i++)
    for(int j = 0; j < (1 << n); j++)
        if(j & (1 << i)){
            f[j] += f[j ^ (1 << i)];//子集和
			//或f[j ^ (1 << i)] += f[j];
        }

时间复杂度$O(n \cdot 2^n)$

求双阶乘的最后64位

跑起来很快, 大概是$O(log^2(n))$级别 \[ (2n - 1)!! = \prod_{i = 1,i += 2}^{2n-1}i \]

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
lld C[64][64], jc[64];

struct poly
{
    lld x[64];

    poly() { for (int i = 0; i < 64; i++) x[i] = 0; }
    poly(lld c1, lld c0) : poly() { x[1] = c1, x[0] = c0; }
    lld operator()(int k) const { return x[k]; }
    lld &operator()(int k) { return x[k]; }

    poly operator*(const poly &b) const
    {
        poly c;
        for (int i = 0; i < 64; i++)
            for (int j = 0; j <= i; j++)
                c(i) += x[j] * b(i-j);
        return c;
    }

    poly moveto(lld n)
    {
        poly ans;
        for (int i = jc[0] = 1; i < 64; i++)
            jc[i] = jc[i-1] * n;
        for (int i = 0; i < 64; i++)
            for (int j = i; j < 64; j++)
                ans(i) += x[j] * jc[j-i] * C[j][i];
        return ans;
    }
};

map<lld, poly> cache;

poly solve(lld n)
{
    if (cache.count(n)) return cache[n];
    poly half = solve(n / 2);
    half = half * half.moveto(n / 2);
    if (n & 1) half = half * poly(2, n*2-1);
    return cache[n] = half;
}

void init()
{
    poly Px(0, 1);
    cache[0] = Px;

    C[0][0] = C[1][0] = C[1][1] = 1;
    for (int n = 2; n < 64; n++)
    {
        C[n][0] = C[n][n] = 1;
        for (int m = 1; m < n; m++)
            C[n][m] = C[n-1][m-1] + C[n-1][m];
    }
}

int main()
{
    init();
    int T; lld n;
    scanf("%d", &T);

    while (T--)
    {
        scanf("%lld", &n);
        printf("%llu\n", solve(n)(0));
    }

    return 0;
}

\[ \sum_{i = 1}^n \sum_{c \in{cnt} \and c < t[i]}\frac{(n - i)!}{(\prod_{c' \in cnt \and c' \neq c}c'!) \cdot (c-1)!} \]

solution-CF1715E

Posted on 2022-08-22 In 题解

题意

做法

斜率优化可以解决如下问题:

给出一组直线$y_i = k_i x_ + b_i$. 对于任意$x$, 求$y_i$的极大/极小值.

[GAMES104]-作业2

Posted on 2022-05-29 In GAMES104

因为一直要赶一堆莫名其妙的ddl和考试 task3没做日后再补_(:з」∠)_

Color Grading与LUT

调色是增强画面表现力的一种很有效的手段. 本次作业我们处理的颜色在sRGB空间内, 红、绿、蓝三个颜色分别用一个$[0,1]$的实数来表示. 将三原色分别用一个坐标轴来表示, 便构成了如图所示的色彩空间(图片来自defold.com):

color cube

在进行调色时, 我们对渲染结果的每个像素进行处理, 根据像素原本的RGB值将其映射到一个新的色彩空间. 这种映射一般采用LUT(look up table)图来进行. LUT有1D和3D之分, Pilot Engine使用的是3D的LUT, 它看上去像这个样子:

这张图片的长为$1024$个像素, 宽为$32$个像素, 可以看出它由$32$个小色块组成. 将这些小方块自上而下拼接在一起, 便组成了一个$32 \times 32 \times 32$的色彩空间. 我们需要做的便是根据每个像素的RGB分量找到它在这张图上对应的像素. 具体来说, 原像素的$G$分量指示了新像素在LUT上的纵坐标, $B$分量指示了新像素所属的色块编号, $R$分量指示了新像素在某个色块内部的横坐标.

编写Shader

我们需要编写的Shader文件是color_grading.frag, 可以在\engine\shader\glsl\中找到.

这个shader已经提供了原像素的颜色in_color和LUT的2D采样器color_grading_lut_texture_sampler. 我们需要根据原像素的颜色计算出调色后的颜色并将其赋值给out_color作为输出.

在计算前, 需要获取像素的颜色值:

1	highp vec4 color = subpassLoad(in_color).rgba;

颜色值是一个四维向量, 四个维度分别代表R、G、B和alpha, 均为$[0,1]$的实数.

在计算出坐标后, 需要从LUT上获取对应的颜色:

1	highp vec4 color_sampled = texture(color_grading_lut_texture_sampler, uv);

其中color_grading_lut_texture_sampler已经提供, uv是一个二维向量, $u$代表横坐标, $v$代表纵坐标, 坐标原点在图片的左上角, $x$轴方向向右, $y$轴方向向下. $u,v$也是$[0,1]$上的实数.

我们还需要获取LUT的长宽:

1	highp ivec2 lut_tex_size = textureSize(color_grading_lut_texture_sampler, 0);

有了如上函数, 不难计算出每个像素在LUT上的对应位置. 设$C$是LUT的宽, $R,G,B$分别为像素的原颜色分量, 则有: \[ u = \frac{\lfloor B \cdot C \rfloor + R}{C}\\ v = G \] 为了提高准确度和防止越界, 可以在LUT的边界上各留出0.5个像素, 将RGB值映射到$C-1$个像素中去. 具体细节见代码:

#version 310 es

#extension GL_GOOGLE_include_directive : enable
//#extension GL_KHR_vulkan_glsl : enable
#include "constants.h"

layout(input_attachment_index = 0, set = 0, binding = 0) uniform highp subpassInput in_color;

layout(set = 0, binding = 1) uniform sampler2D color_grading_lut_texture_sampler;

layout(location = 0) out highp vec4 out_color;
void main()
{
    highp ivec2 lut_tex_size = textureSize(color_grading_lut_texture_sampler, 0);
    highp float COLORS       = float(lut_tex_size.y);//LUT的每一维的颜色数量
    highp vec4  color        = subpassLoad(in_color).rgba;

    highp float max_color = COLORS - 1.0;//允许使用的最大元素
    highp float half_pix_u = 0.5 / float(lut_tex_size.x);//预留出半个像素
    highp float half_pix_v = 0.5 / float(lut_tex_size.y);

    highp float part = floor(color.b * max_color);//根据B计算出该像素所属的色块编号
    highp float threshold = max_color / COLORS;//现在只有COLORS-1个颜色可用
    highp float u = (half_pix_u + threshold * color.r + part) / COLORS;
    highp float v = (half_pix_v + threshold * color.g);
    //预留出0.5像素的边界.例如假设COLORS = 32,在G = 0时v = 0.5/32,G = 1时v = 31.5/32

    highp vec2 uv = vec2(u,v);
    
    highp vec4 color_sampled = texture(color_grading_lut_texture_sampler, uv);

    out_color =  color_sampled;
}

效果展示与自定义LUT

在\engine\asset\texture\lut中存放着不同的LUT. 通过修改\engine\asset\global\rendering.global.json中"color_grading_map"项可以替换不同的LUT.

在不进行调色时的效果是这样的:

调色后的效果:

要创建自定义的LUT, 只需在GIMP或者PhotoShop等工具编辑原始LUT图像(可以附一张截图方便观察调色效果).

负片效果的LUT:

偏冷色调(好像太绿了):

PS:调色后会出现看上去很奇怪的色块, 目前还不知道如何解决

[笔记]从0开始搓一个简单的render

Posted on 2022-04-30 Edited on 2022-06-30 In 图形学

本文的主要内容是对Dmitry V. Sokolov的tinyrender课程的记录和补充.

Step 0：环境搭建

我们需要tgaimage.h来操作TGA文件, model.h来操作模型. 头文件和实现可以在这里找到:

https://github.com/ssloy/tinyrenderer/tree/f6fecb7ad493264ecd15e230411bfb1cca539a12

可以用如下代码来测试:

#include "tgaimage.h"
const TGAColor white = TGAColor(255, 255, 255, 255);
const TGAColor red   = TGAColor(255, 0,   0,   255);
int main(int argc, char** argv) {
        TGAImage image(100, 100, TGAImage::RGB);
        image.set(52, 41, red);
        image.flip_vertically(); // i want to have the origin at the left bottom corner of the image
        image.write_tga_file("output.tga");
        return 0;
}

得到一个有红点的TGA文件:

Step 1：画一条直线

朴素算法

我们需要在二维位图上画一条从$(x_0,y_0)$到$(x_1,y_1)$的直线.

直线的两点式为: \[ \frac{y - y_0}{y_1 - y_0} = \frac{x - x_0}{x_1 - x_0} \] 因而对任一点$x$, 有 \[ y = \frac{x - x_0}{x_1 - x_0} \cdot(y_1 - y_0) + y_0 \] 我们只要枚举$x$, 计算出每个$x$对应的$y$即可画出直线. 需要注意的是, 这个算法默认$(x_1,y_1)$在$(x_0,y_0)$的右上方, 并且绘制时需要枚举$\Delta x,\Delta y$中较大的一方以保证线段不会断开.

代码如下:

void Drawer::Line(int x0, int y0, int x1, int y1, TGAColor color) {
    bool flip = 0;
    if(abs(x0 - x1) < abs(y0 - y1)){
        flip = 1;
        std::swap(x0,y0),std::swap(x1,y1);
    }
    if(x0 > x1){
        std::swap(x0,x1),std::swap(y0,y1);
    }
    for(int x = x0;x <= x1;x++){
        double t = (double)(x - x0)/(x1 - x0);
        int y = t * (y1 - y0) + y0;
        if(!flip)set(x,y,color);
        else set(y,x,color);
    }
}

Bresenham直线算法

朴素算法已经能够满足绘图需要, 但是使用了浮点乘法除法, 效率很低.

我们设直线的斜率$k = \frac{y_1 - y_0}{x_1 - x_0}$, 有 \[ y = (x - x_0) \cdot k + y_0 \] 我们假设当前画笔所在位置的是$(x,y_c)$, 误差为$E = y- y_c$. 当$x$加$1$时, $E$增加$k$. 当$E \ge \frac{1}{2}$时, 我们需要将$y_c$加或减$1$同时$E$减去$1$. 这样便避免了每次都要进行的浮点除法运算.

具体来说, 我们只需将不等式的两边乘上$2 \cdot (x_1 - x_0)$即可. 此时误差每次增加$2 \cdot (y_1 - y_0)$. 当$E \ge (x_1 - x_0)$时需要将$y_c$加或减$1$同时$E$减去$2 \cdot (x_1 - x_0)$.

void Drawer::Line(int x0, int y0, int x1, int y1, TGAColor color) {
    bool flip = 0;
    if(abs(x0 - x1) < abs(y0 - y1)){
        flip = 1;
        std::swap(x0,y0),std::swap(x1,y1);
    }
    if(x0 > x1){
        std::swap(x0,x1),std::swap(y0,y1);
    }
    int dx = x1 - x0,dy = y1 - y0;
    int e = 0,de = 2*abs(dy);
    int y = y0;
    for(int x = x0;x <= x1;x++){
        e += de;
        if(e >= dx) {
            e -= dx*2;
            y += (dy > 0?1:-1);
        }
        if(!flip)set(x,y,color);
        else set(y,x,color);
    }
}

Step 2：画一个三角形

扫描线法

我们先用一种十分古老的方法来画一个三角形: 用水平或者竖直的线一行一行地绘制.

在绘制前, 先对三角形的三个顶点按横坐标递增进行排序, 然后标号为$P_0,P_1,P_2$(以下默认左下角为坐标原点). 同时将边$P_0P_1$标为$L_0$, $P_1P_2$标为$L_1$, $P_2P_3$标为$L_2$, 如图:

标号后的三角形可按$x = x_1$这条直线分割成两部分. 对两部分分别进行绘制即可. 注意处理一下斜率不存在时的细节.

void Drawer::Triangle(Vec2i t0, Vec2i t1, Vec2i t2, TGAColor color) {
    int x0,y0,x1,y1,x2,y2;
    std::vector<std::pair<int,int> >top = {{t0.u,t0.v},{t1.u,t1.v},{t2.u,t2.v}};
    std::sort(top.begin(),top.end());
    std::tie(x0,y0) = top[0];
    std::tie(x1,y1) = top[1];
    std::tie(x2,y2) = top[2];
    int mx = std::max(x1,x2);
    double k0 = (double)(y1 - y0)/(x1 == x0?1:x1 - x0);
    double k1 = (double)(y2 - y1)/(x2 == x1?1:x2 - x1);
    double k2 = (double)(y2 - y0)/(x2 == x0?1:x2 - x0);
    for(int x = x0;x <= x2;x++){
        int yc2 = (x - x0) * k2 + y0;
        if(x <= x1){
            int yc0 = (x - x0) * k0 + y0;
            Line(x,yc0,x,yc2,color);
        }
        else{
            int yc1 = (x - x1) * k1 + y1;
            Line(x,yc1,x,yc2,color);
        }
    }
}

重心坐标系

一般定义

设$V_1,...,V_n$是$n$维向量空间$V$中单纯形的顶点, 对于$V$中任意一点$P$, 有 \[ (\sum_{i = 1}^n \lambda_i)\cdot P = \sum_{i = 1}^n \lambda_i \cdot V_i \] 则系数$\lambda_1,...,\lambda_n$称为$P$关于$V_1,...V_n$的重心坐标. 一般规定$\sum_{i = 1}^n \lambda_i = 1$, 此时称为正规化的重心坐标. 以下讨论的重心坐标均为正规化的.

二维平面中的重心坐标系

定义

二维平面上, 单纯形为三角形. 假设给出一三角形$\Delta ABC$和三角形所在平面上的一点$P$, 则$P$可被如下重心坐标表示: \[ P = \alpha A + \beta B + \gamma C \] 考虑下图(来源见链接):

将$\vec{AB},\vec{AC}$作为坐标系, 可以得到 \[ \begin{align} P =& A + \beta(\vec{AB}) + \gamma(\vec{AC})\\ =& A + \beta(B - A) + \gamma(C - A)\\ =& (1 - \beta - \gamma)A + \beta B + \gamma C \end{align} \] 这样便可以从坐标系的角度来解释重心坐标, 同时能得到 \[ \alpha = 1 - \beta - \gamma \]

求解:解方程组

将点带入重心坐标的定义, 可以得到: \[ \left \{ \begin{array}{cccc} x_p = (1 - \beta - \gamma) x_a + \beta x_b + \gamma x_c \\ y_p = (1 - \beta - \gamma) y_a + \beta y_b + \gamma y_c \end{array} \right. \] 提取出$\beta \space \gamma$: \[ \left \{ \begin{array}{cccc} x_p - x_a = \beta(x_b - x_a) + \gamma (x_c - x_a)\\ y_p - y_a = \beta(y_b - y_a) + \gamma (y_c - y_a) \end{array} \right. \] 转化为矩阵形式: \[ \left[ \begin{array}{} x_b - x_a & x_c - x_a\\ y_b - y_a & y_c - y_a \end{array} \right] \left[ \begin{array}{} \beta \\ \gamma \end{array} \right] = \left[ \begin{array}{} x_p - x_a\\ y_p - y_a \end{array} \right] \] 这样便可求解出$\beta \space \gamma$进而求出$\alpha$

求解:向量叉积

把坐标系的表达式变换一下: \[ \beta \vec{AB} + \gamma \vec{AC} + \vec{PA} = 0 \] 拆开后: \[ \left \{ \begin{array}{cccc} \beta\vec{AB}_x + \gamma \vec{AC}_x + \vec{PA}_x = 0\\ \beta\vec{AB}_y + \gamma \vec{AC}_y + \vec{PA}_y = 0 \end{array} \right. \] 转化为矩阵形式: \[ \left[ \begin{array}{} \beta & \gamma & 1 \end{array} \right] \left[ \begin{array}{} \vec{AB}_x\\ \vec{AC}_x\\ \vec{PA}_x \end{array} \right] = 0\\ \left[ \begin{array}{} \beta & \gamma & 1 \end{array} \right] \left[ \begin{array}{} \vec{AB}_y\\ \vec{AC}_y\\ \vec{PA}_y \end{array} \right] = 0\\ \] 从几何意义考虑, 这相当于$[\beta ,\gamma,1]$这个向量和后两个向量分别垂直. 因而我们求出后两个向量的叉积(假设结果为$[x,y,z]$), 则$[\frac{x}{z}, \frac{y}{z},1]$等于$[\beta,\gamma,1]$.

需要注意$z = 0$的特殊情况, 根据叉乘的定义, 此时有$\vec{AB}_x \cdot \vec{AC}_y - \vec{AB}_y \cdot \vec{AC}_x = 0$, 即$\vec{AB} \cdot \vec{AC} = 0$, 说明$ABC$三点共线. 这时返回一个任意的负重心坐标即可.

求解:三角形面积

如图, 假设三角形面积分别为$S_A,S_B,S_C$. 则 \[ \alpha = \frac{S_A}{S_A+S_B+S_C}\\ \beta = \frac{S_B}{S_A+S_B+S_C}\\ \gamma = \frac{S_C}{S_A+S_B+S_C} \]

应用

重心坐标有很多应用, 在这里我们能用到的结论是: 若$\alpha,\beta,\gamma$均大于$0$, 则点$P$在三角形$\Delta ABC$的内部; 若$\alpha,\beta,\gamma$有一个等于$0$, 则点$P$在三角形的边上; 若$\alpha,\beta,\gamma$有两个等于$0$, 则点$P$在三角形的顶点上.

在绘图时, 用一个最小的矩形包围住三角形, 再枚举并判断每个点是否在三角形内部即可. 下面的代码采用了向量叉积法:

Vec3f Drawer::Barycentric(const Vec2i (&vertex)[3],const Vec2i & P){
    Vec3f result =  //cal cross
            Vec3f(vertex[1].x - vertex[0].x,vertex[2].x - vertex[0].x,vertex[0].x - P.x) ^
            Vec3f(vertex[1].y - vertex[0].y,vertex[2].y - vertex[0].y,vertex[0].y - P.y);
    if(fabs(result.z) < 1)return Vec3f(-1,1,1);
    result.x /= result.z,result.y /= result.z;
    return {(float)1.0 - result.x - result.y,result.x,result.y};
}
void Drawer::Triangle(const Vec2i (&vertex)[3],const TGAColor & color) {
    int lx = std::max(std::min({vertex[0].x,vertex[1].x,vertex[2].x}),0),rx = std::min(std::max({vertex[0].x,vertex[1].x,vertex[2].x}),get_width()-1);
    int ly = std::max(std::min({vertex[0].y,vertex[1].y,vertex[2].y}),0),ry = std::min(std::max({vertex[0].y,vertex[1].y,vertex[2].y,}),get_height()-1);
    for(int x = lx;x <= rx;x++){
        for(int y = ly;y <= ry;y++){
            Vec3f bary = Barycentric(vertex,Vec2i(x,y));
            bool out = (bary.x < 0 or bary.y < 0 or bary.z < 0);
            if(!out)set(x,y,color);
        }
    }
}

使用tinyrender提供的模型来测试一下我们的三角形绘制程序(先给每个面一个随机颜色):

加上一点阴影并剔除背面的三角形:

Step 3 ：z-buffer

现在我们需要考虑剔除掉不能被我们看见的像素(Step 2通过Back-face culling剔除了一些面, 但是仍然不正确, 例如嘴部整个消失了). z-buffer的思想是维护一个像素点在待渲染平面上的一个距离buffer, 每次用离相机更近的点来剔除已经存在的点. z-buffer的实现非常简单, 下面主要讨论如何求每个点距相机的距离.

首先考虑一维的情形. 假设我们要绘制从$A(x_0,y_0)$到$B(x_1,y_1)$的一条线段, 当前绘制到点$P(x,y)$, 摄像机的视角与$y$轴的反方向相同:

在绘制时我们要通过$x$坐标计算出$y$来. 不难得到 \[ \frac{y - y_0}{y1 - y_0} = \frac{x - x_0}{x_1 - x_0} \] 我们设$\beta = \frac{x - x_0}{x_1 - x_0}$, 则有 \[ y = (1 - \beta)y_0 + \beta y_1 \] 是不是很熟悉? $\alpha = 1 - \beta$和$\beta$就是$P$关于线段$AB$的重心坐标(一维单纯形是线段). 我们将其推广到二维, 便能求出三角形每个点距离相机的距离$z$: \[ z = \alpha z_0+ \beta z_1 + \gamma z_2 \]

void Drawer::Triangle(const Vec3i (&vertex)[3],const TGAColor & color,double * zbuffer) {
    const double eps = -1e-6;//防止三角形接合处有未绘制的点
    int lx = std::max(std::min({vertex[0].x,vertex[1].x,vertex[2].x}),0),rx = std::min(std::max({vertex[0].x,vertex[1].x,vertex[2].x}),get_width()-1);
    int ly = std::max(std::min({vertex[0].y,vertex[1].y,vertex[2].y}),0),ry = std::min(std::max({vertex[0].y,vertex[1].y,vertex[2].y,}),get_height()-1);
    for(int x = lx;x <= rx;x++){
        for(int y = ly;y <= ry;y++){
            Vec3f bary = Barycentric(vertex,Vec2i(x,y));
            if(bary.x < eps or bary.y < eps or bary.z < eps)continue;
            int idx = x + y * get_width();
            double z = bary.x * vertex[0].z + bary.y * vertex[1].z + bary.z * vertex[2].z;
            if(zbuffer[idx] < z){
                zbuffer[idx] = z;
                set(x,y,color);
            }
        }
    }
}

效果如下:

Step 4: 透视投影(WIP)

这一部分tinyrender讲得很含糊而且不直观, 所以直接按照GAMES101的内容去实现. 本节只讲述实现. obj文件格式可以在https://en.wikipedia.org/wiki/Wavefront_.obj_file找到.

透视投影的实现似乎有问题, 先鸽了

Step 5:相机

写的不是很满意, 也鸽了

[题解]CF1642F - Two Arrays

Posted on 2022-03-18 In 题解

感觉这场F比E简单多了= =

题意

给出一个$n \times m$的数组 $(1 \leq n \leq 10^5,1 \leq m \leq 5)$. 数组的每一行元素均互不相同, 且每一行有一个权值$w$ $(1 \leq w \leq 10^9)$.

你需要找出不同的两行, 使得这两行在不存在重复元素的同时权值最小. 输出权值, 无解输出-1.

题解

暴力做需要$O(n^2m)$.

将行按权值排序, 把每一个元素对应的行用bitset存下来, 0代表非法1代表合法. 再在枚举每行的时候找到第一个不冲突的行更新答案就可以做到$O(\frac{n^2m}{w})$, 其中$w = 32$或$64$. 但是这样需要$O(\frac{n^2m}{8})$的空间, 会MLE.

考虑分块. 设块大小为$B$, 将每个元素按出现次数分类:

出现次数大于$B$, 这部分元素用bitset存下来, 在枚举每行时使用按位与来更新, 时间复杂度$O(\frac{n^2m}{Bw})$, 空间复杂度$O(\frac{n^2m}{8B})$.
出现次数小于$B$, 这部分元素在枚举每行时直接更新, 时间复杂度$O(B^2)$, 空间复杂度$O(n)$.

(感觉复杂度算的很不对劲啊- - 反正能过)

然后B选一个差不多的大小比如$\sqrt{nm}$就能过了.

代码

//
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define endl '\n'
#define ALL(x) (x).begin(),(x).end()
#define PRINT(___) for(auto ____:(___)){cout << ____ << " ";}cout << endl;
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
ll Pow(ll x,ll d,ll p){ll ta=1;if(d==0)return 1%p;x %= p;ll a=Pow(x,d/2,p);ta=a*a%p;if(d%2)ta=ta*x%p;return ta%p;}
const int maxn = 1e5+10;//CHANGE IT!
const int maxm = maxn*2;
const int p = 1e9+7;//CHANGE IT!
int raw[maxn],mp[maxn * 5];
const int maxb = 708;
const int maxc = 100006;
bitset<maxc>b[maxb + 10];
vector<int>ban[maxn * 5];
struct Node{
	vector<int>r;
	int w;
	Node(vector<int> r = {},int w = 0):r(r),w(w){}
	bool operator < (const Node & x)const{
		return w < x.w;
	}
}node[maxn];
int cnt[maxn * 5];
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.setf(ios::fixed,ios::floatfield);
	cout.precision(12);	 
	int n,m;
	cin >> n >> m;
	map<int,int>hsh;
	for(int i = 1;i <= n;i++){
		vector<int>r(m);
		int w;
		for(int j = 0;j < m;j++){
			cin >> r[j];
			hsh[r[j]]++;
		}
		cin >> w;
		node[i] = Node(r,w);
	}
	int tmp = 1;
	for(auto & it:hsh){
		cnt[tmp] = it.second;
		it.second = tmp++;
	}
	sort(node+1,node+n+1);
	int ans = 2e9 + 10,idx = 1;
	for(int i = 1;i <= n;i++){
		for(int j = 0;j < m;j++){
			int & v = node[i].r[j];
			v = hsh[v];
			if(cnt[v] > maxb){
				if(!mp[v]){
					b[idx].set();
					mp[v] = idx++;
				}
				b[mp[v]].reset(i);
			}
			else ban[v].push_back(i);
		}
	}
	for(int i = 1;i <= n;i++){
		bitset<maxc>vai;
		vai.set();
		vai.reset(0);
		for(int j = 0;j < m;j++){
			int v = node[i].r[j];
			if(cnt[v] > maxb){
				vai &= b[mp[v]];
			}
			else{
				for(int x:ban[v])vai.reset(x);
			}
		}
		int x = vai._Find_first();
		if(x <= n)ans = min(ans,node[i].w + node[x].w);
	}
	if(ans == 2e9+10)ans = -1;
	cout << ans;
	return 0;
}

[题解]Gym103495-H Reverse the String

Posted on 2022-03-11 Edited on 2022-03-18 In 题解

中规中矩的字符串题.

题意

给出一个字符串$s$, 你可以选取$s$的任意子串将其翻转得到$s'$, 该操作最多进行一次. 求字典序最小的$s'$.

$|s| \leq 10^5, \sum |s| \leq 1.5 \cdot 10 ^ 6$

题解

以下字符串的下标均从$0$开始. $lcp(s,t)$代表$s$和$t$的最长公共前缀. $rev(s)$表示将$s$翻转后的字符串.

part 1：优化做法

设$n = |s|$, 朴素的暴力需要至少$O(n^2)$来枚举翻转子串的两个端点$L,R$. 我们需要想办法固定住其中的一个端点.

首先尝试放开限制, 把翻转变成任意排列. 这样将$s$排序后得到的$s'$一定是最优的. 记这样的$s'$为$t$.

回到原题, 容易看出, $s$和$t$的最长公共前缀部分无论我们如何操作都不可能再让字典序变小, 而在这之后的部分是有可能减小字典序的. 且因为字典序的性质(比较第一个不相同的字符), 我们必须尝试让$s'$和$t$不同的第一个字符最小. 因此$L = lcp(s,t)$.

part 2：快速比较两字符串大小

设当前的最优答案为$ans$, 在固定了$L$之后, 只需$O(n)$枚举$R$, 并快速比较$ans$和$s'$的大小即可. 这里可以用$lcp$来比较两个字符串的大小.

设需要比较$A = s[a...b]$和$B = s[c...d]$的关系. 若$lcp(a,c) \ge \min(|A|,|B|)$,则$A < B$等价于$|A| < |B|$ 否则, $A < B$等价于$rank[a] < rank[c]$

int CompString(int s1,int len1,int s2,int len2){
    	//s:子串起始位置 len:子串长度
		int tlcp = (s1 == s2?min(len1,len2):LCP(s1,s2));
		if(tlcp >= min(len1,len2)){
			return len1 < len2?-1:(len1 == len2?0:1);
		}
		else return c[s1] < c[s2]?-1:(c[s1] == c[s2]?0:1);
	}

注意: 可以只用$lcp$来比较字符串大小(这样就可以用hash等方法来求$lcp$而不用写SA). 只要把上面的第6行替换为

1	else return s[s1 + tlcp] < s[s2 + tlcp]?-1:(s[s1 + tlcp] == s[s2 + tlcp]?0:1);

即可.

part 3 ：细节部分

首先将整个翻转后的$s$拼接在原字符串后, 中间用特殊字符隔开, 例如u = s + "$" + rev(s). 再使用后缀数组处理$u$. 这样在原串中位置为$i$的字符, 在反串中的位置就是$2*n - i$.

假设当前枚举到$R$, 最优答案为$R'$, 如下图:

我们需要比较的是$s[L,R'-1] + rev(s[R',R])$和$rev(s[L,R])$的大小, 如果后者更小则更新答案. 为了方便比较, 可以将每段字符划分成AB两段:

现在需要比较的就是$rev(A_1) + B_1$和$rev(A_2) + rev(B_2)$的大小. 由于A, B长度分别相等, 可以用part2的方法分别比较$rev(A_1)$与$rev(A_2)$, $B_1$与$rev(B_2)$的大小

每段的(起始位置,长度)可以通过简单推导得到:

$rev(A_1)$: $(2*n - R',R' - L + 1)$

$rev(A2)$: $(2*n - R,R - L + 1)$

$B_1$: $(R' + 1,i - R')$

$rev(B_2)$: $(2*n - R + R' - L + 1,i - R')$

其他细节见代码. 时间复杂度$O(n \log n)$. jls好像有$O(n)$做法, 有没有好哥哥教教QAQ

代码

//
#pragma GCC optimize("Ofast,unroll-loops")
#pragma GCC target("avx,avx2,sse,sse2")
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define endl '\n'
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
const int maxn = 2e5+10;//CHANGE IT!
const int maxm = maxn*2;
const int p = 1e9+7;//CHANGE IT!
ll Pow(ll x,ll d,ll p){ll ta=1;if(d==0)return 1%p;x %= p;ll a=Pow(x,d/2,p);ta=a*a%p;if(d%2)ta=ta*x%p;return ta%p;}
struct Node{
	int a,b;
	Node(int a = 0,int b = 0):a(a),b(b){}
	bool operator < (const Node & x){
		return a == x.a?b < x.b:a < x.a;
	}
}; 
int lcp[maxn],ST[maxn][20];
struct SA{
	string s;
	int n;
	vector<int>sa,c;
	SA(string s = ""):s(s){
		sa.clear(),c.clear();
		s += '$';
		n = s.size();
		sa.assign(n,0);
		c.assign(n,0); 
	}
	void count_sort(vector<int> & sa,vector<int> & c){
		vector<int>pos(n),cnt(n),sa_new(n);
		for(int x:c)cnt[x]++;
		for(int i = 1;i < n;i++)pos[i] = pos[i-1] + cnt[i-1];
		for(int x:sa){
			int d = c[x];
			sa_new[pos[d]] = x;
			pos[d]++;
		}
		sa = sa_new;
	}
	void cal(){
		//init
		{
			vector<Node>A(n);
			for(int i = 0;i < n;i++)A[i] = Node(s[i],i);
			sort(A.begin(),A.end());
			for(int i = 0;i < n;i++)sa[i] = A[i].b;
			c[sa[0]] = 0;
			for(int i = 1;i < n;i++){
				c[sa[i]] = c[sa[i-1]] + (A[i].a != A[i-1].a);
			}
		}
		int k = 0;
		while((1<<k) < n){
			for(int i = 0;i < n;i++)sa[i] = (sa[i] - (1<<k) + n)%n;
			count_sort(sa,c);
			vector<int>c_new(n);
			c_new[sa[0]] = 0;
			for(int i = 1;i < n;i++){
				pair<int,int>now = {c[sa[i]],c[(sa[i] + (1<<k)) % n]};
				pair<int,int>prev = {c[sa[i-1]],c[(sa[i-1] + (1<<k)) % n]};
				c_new[sa[i]] = c_new[sa[i-1]] + (now != prev);
			}
			k++;
			c = c_new;
		}
	}
	void ST_init(){
		for(int i = 1;i <= n;i++)ST[i][0] = lcp[i];
		for(int j = 1;(1<<j) <= n;j++)
			for(int i = 1;i + (1<<j) -1 <= n;i++)
				ST[i][j] = min(ST[i][j - 1],ST[i + (1<<(j - 1))][j - 1]);
	}
	int LCP(int L,int R){
		L = c[L],R = c[R];
		if(L > R)swap(L,R);
		L++;
		int k = log2(R - L + 1);
		return min(ST[L][k],ST[R - (1<<k) + 1][k]);
	}
	void getlcp(){
		int k = 0;
		for(int i = 0;i < n-1;i++){
			int j = sa[c[i] - 1];
			while(s[i+k] == s[j+k])k++;
			lcp[c[i]] = k;
			if(k)k--;
		}
		ST_init();
	} 
	int CompString(int s1,int len1,int s2,int len2){
		int tlcp = (s1 == s2?min(len1,len2):LCP(s1,s2));
		if(tlcp >= min(len1,len2)){
			return len1 < len2?-1:(len1 == len2?0:1);
		}
		else return s[s1 + tlcp] < s[s2 + tlcp]?-1:(s[s1 + tlcp] == s[s2 + tlcp]?0:1);
	}
}sa;
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.setf(ios::fixed,ios::floatfield);
	cout.precision(12);	 
	int __;
	cin >> __;
	while(__--){
		string s;
		cin >> s;
		string u = s;
		sort(u.begin(),u.end());
		int L = -1;
		int n = s.size();
		for(int i = 0;i < n;i++){
			if(s[i] != u[i]){
				L = i;
				break;
			}
		}
		if(L == -1){
			cout << s << endl;
			continue;
		}
		u = s;
		reverse(u.begin(),u.end());
		s = s + "$" + u;
		sa = SA(s);
		sa.cal();
		sa.getlcp();
		int R = L;
		for(int i = L;i < n;i++){
			if(sa.CompString(L,i - L + 1,2*n - i,i - L + 1) >= 0){
				int x = sa.CompString(2*n - R,R - L + 1,2*n - i,R - L + 1);
				if(x > 0){
					R = i;					
				}
				else if(x == 0 and sa.CompString(R + 1,i - R,2*n - i + R - L + 1,i - R) >= 0){
					R = i;
				}
			}
		}
		reverse(u.begin(),u.end());
		reverse(u.begin() + L,u.begin() + R + 1);
		cout << u << endl;
	}
	return 0;
}

solution-CF1625E1

Posted on 2022-01-22