In mathematics you don't understand things, you just get used to them.

分类 笔记 下的文章

区间加法可以差分的原因是 $\Delta_i = \Delta_{i-1}$,所以令 $\textit{d}_i = \Delta_i-\Delta_{i-1}$,对于一个一般点的 $\Delta_i = f(\{\Delta_1, \dots, \Delta_{i-1}\})$,我们可以令 $\textit{d}_i = \Delta_i-f(\{\Delta_1, \dots, \Delta_{i-1}\})$。

这得出了广义一点的差分。那么前缀和的方式即 $\textit{d}_i \gets d_i+f(\{d, \dots, d_{i-1}\})$(其实这里的 $f$ 改为下标的函数可能更准确一点)。

举个例子,对区间 $[l, r]$ 中的 $i \in [l, r]$ 操作 $a_i \gets a_i+f_{i-l+1}$,其中 $f$ 是 Fibonacci 数列,也就是说 $f_i = f_{i-1}+f_{i-2}$,差分数组即 $d_i = a_i-a_{i-1}-a_{i-2}$,操作即 $d_l \gets d_l+1$,$d_{r+1} \gets d_{r+1}+f_{r-l+2}$,$d_{r+2} \gets d_{r+2}+f_{r-l+1}$,前缀和方法即 $d_i \gets d_i+d_{i-1}+d_{i-2}$。

差分是导函数的一个演绎,即令研究范围中的 $\epsilon = 1$ 而得来的(在实数中令 $\epsilon = \lim\limits_{\frac{1}{x} \rightarrow +\infty}x$ 做差分即可得到导函数)。

Quack 知道好多东西,把它们都做成 ppt。

inv_gcd 还可以用递推矩阵算。

void exgcd(int a, int b){
  r[0] = a, r[1] = b;
  int i = 2;
  for ( ; r[i - 1]; i++) {
    r[i] = r[i - 2] % r[i - 1];
    q[i - 1] = r[i - 2] / r[i - 1];
  }
  d = r[i -= 2];
  mat res(1, 0, 0, 1);
  for (int j = i - 1; j; --j)
    res.mul(mat(0, 1, 1, -q[j]));
  s = res.get(1, 0); t = res.get(1, 1);
}

还有个 trick 叫 binary gcd。

int gcd(int a, int b){
  if (a == b) return a;
  if ((a & 1) && (b & 1)) return gcd(b, abs(a - b));
  else if ((a & 1) && (!(b & 1))) return gcd(a, b >> 1);
  else if ((!(a & 1)) && (b & 1)) return gcd(a >> 1, b);
  else return gcd(a >> 1, b >> 1) << 1;
}

复杂度和正确性平凡吧,在写高精度 gcd 的时候可以考虑。接下来看一个典中典问题:

证明:$(a^n-1, a^m-1) = a^{(n, m)}-1$。

P.f.:不妨设 $n \leqslant m$,首先 $(a^n-1, a^m-1) = (a^n(a^{m-n}-1), a^n-1) = (a^{m-n}-1, a^n-1)$,一直迭代可以得到 $(a^n-1, a^m-1) = (a^{(n, m)}-1, a^{(n, m)}-1) = a^{(n, m)}-1$。

然后这个有个 more general 的版本,就是对于任意整数 $a, b$,满足 $(a, b) = 1$ 则有 $(a^n-b^n, a^m-b^m) = a^{(n, m)}-b^{(n, m)}$,证明是类似的。

但是如果说我们还有一个 more more general 的版本呢?

Theorem:若 $f_n$ 是一个整数序列,并满足 $f_0 = 0$,$f_{n} \equiv f_{n-m} \pmod {f_m},n > m$,有 $(f_n, f_m) = f_{(n, m)}$。

P.f.:考虑 induction,当 $n = m$,$n = 0$,$m = 0$ 时,上述命题为一个,一个一个一个,一个一个一个一个一个一个一个

什么栈什么溢什么出上还有一些 q-analog 的高论,不过我连上面的证明都还没看懂。引一下原 writer:

Mimic in expts a subtractive Euclidean algorithm $\rm\,(n,m) = (\color{#0a0}{n\!-\!m},m)$

$$\begin{align} \rm{e.g.}\ \ &\rm (f_5,f_2) = (f_3,f_2) = (f_1,f_2) = (f_1,f_1) = (f_1,\color{darkorange}{f_0})= f_1 = f_{\:\!(5,\,2)}\[.3em]
{\rm like}\ \ \ &(5,\ 2)\, =\:\! (3,\ 2)\, =\:\! (1,\ 2)\:\! =\:\! (1,\ 1)\:\! =\:\! (1,\ \color{darkorange}0)\:\! = 1,\ \ {\rm since}\end{align}\qquad$$

$\rm\ f_{\,n}\: :=\ a^n\!-\!1\ =\ a^{n-m} \: \color{#c00}{(a^m\!-\!1)} + \color{#0a0}{a^{n-m}\!-\!1},\,\ $ hence $\rm\:\ {f_{\,n}\! = \color{#0a0}{f_{\,n-m}}\! + k\ \color{#c00}{f_{\,m}}},\,\ k\in\mathbb Z,\:$ thus

Theorem $\: $ If $\rm\ f_{\, n}\: $ is an integer sequence with $\rm\ \color{darkorange}{f_{0} =\, 0},\: $ $\rm \:{ f_{\,n}\!\equiv \color{#0a0}{f_{\,n-m}}\ (mod\ \color{#c00}{f_{\,m})}}\ $ for all $\rm\: n > m,\ $ then $\rm\: (f_{\,n},f_{\,m})\ =\ f_{\,(n,\:m)}, \: $ where $\rm\ (i,\:j)\ $ denotes $\rm\ gcd(i,\:j).\:$

Proof $\ $ By induction on $\rm\:n + m\:$. The theorem is trivially true if $\rm\ n = m\ $ or $\rm\ n = \color{darkorange}0\ $ or $\rm\: m = \color{darkorange}0.\:$
So we may assume $\rm\:n > m > 0\:$.$\ $ Note $\rm\ (f_{\,n},f_{\,m}) = (\color{#0a0}{f_{\,n-m}},\color{#c00}{f_{\,m}})\ $ follows by $\rm\color{#90f}{Euclid}$ & hypothesis.
Since $\rm\ (n-m)+m \ <\ n+m,\ $ induction yields $\rm\, \ (f_{\,n-m},f_{\,m})\, =\, f_{\,(n-m,\:m)} =\, f_{\,(n,\:m)}.$

$\rm\color{#90f}{Euclid}\!:\ A\equiv a\pmod{\! m}\,\Rightarrow\ (A,m) = (a,m)\,$ is the reduction (descent) step used both above and in the Euclidean algorithm $\rm\: (A,m) = (A\bmod m,\,m),\, $ the special case $\,\rm f_{\:\!n} = n\,$ above.

This is a prototypical strong divisibility sequence. Same for Fibonacci numbers.


Alternatively it has a natural proof via the Order Theorem $\ a^k\equiv 1\iff {\rm ord}(a)\mid k,\,$ viz.

$$\begin{eqnarray}\ {\rm mod}\:\ d\!:\ a^M\!\equiv 1\equiv a^N&\!\iff\!& {\rm ord}(a)\ |\ M,N\!\color{#c00}\iff\! {\rm ord}(a)\ |\ (M,N)\!\iff\! \color{#0a0}{a^{(M,N)}\!\equiv 1}\[.2em]
{\rm i.e.}\ \ \ d\ |\ a^M\!-\!1,\:a^N\!-\!1\! &\!\iff\!\!&\ d\ |\ \color{#0a0}{a^{(M,N)}\!-\!1},\qquad\,\ {\rm where} \quad\! (M,N)\, :=\, \gcd(M,N)
\end{eqnarray}\ \ \ \ \ $$

Thus, by above $\, a^M\!-\!1,\:a^N\!-\!1\ $ and $\, a^{(M,N)}\!-\!1\ $ have the same set $\,S\,$ of common divisors $\,d,\, $ therefore they have the same greatest common divisor $\ (= \max\ S).$

Note We used the GCD universal property $\ a\mid b,c \color{#c00}\iff a\mid (b,c)\ $ [which is the definition of a gcd in more general rings]. Compare that with $\ a<b,c \!\iff\! a< \min(b,c),\, $ and, analogously, $\,\ a\subset b,c\iff a\subset b\cap c.\ $ Such universal "iff" characterizations enable quick and easy simultaneous proof of both directions.

The conceptual structure that lies at the heart of this simple proof is the ubiquitous order ideal. $\ $ See this answer for more on this and the more familiar additive form of a denominator ideal.

link。

一个 dp(拜谢 ly)和切入点都略有不同的做法,并不需要观察啥性质。

原问题针对子序列进行规划,自然地想到转而对前缀进行规划。接下来我们考虑一个前缀 $[1, i]$ 以及一个 $j \in [1, i]$ 对答案的贡献,可以写出 $\displaystyle \textit{cont}(j \text{ to } [1, i]) = [\max_{i < k} a_k > a_j] \times \text{the number of LISs containing } j \text{ indexed in } [1, i]$。

这个可以做两个 dp 解决,第一个好解决的静态 dp,即结束在 $j$ 的 LIS 方案数;第二个 dp 有些烦:$j$ 在动,我们考虑的前缀 $[1, i]$ 也在移动。

到这里其实换一下着手点第二个 dp 就变成静态的了,即考虑位置 $j$,直接算 $(j, i)$ 的贡献即可,其中 $i$ 是「满足 $a_i > a_j$ 的最大的 $i$」。于是我们的第二个 dp 就可以被描述为从 $j$ 开始,结束在 $i$ 之前(不包括,因为要保证 $\max_{i < k} a_k > a_j$)的 LIS 方案数。答案即 $\displaystyle \sum_{i=1}^n dp_i \times dp'_i$。

第二个 dp 的具体实现,还需要一个辅助的 dp,其定义是第二个 dp 的定义中去掉贡献范围上界(即 $i$),转移画画图就能理解了。

用下 Fenwick Tree 啥的就能 $O(n \log_2 n)$ 了。

using modint = modint1000000007;
int n, a[200100], pos[200100], id[200100];
modint dp[200100], dp2[200100], dp3[200100], bit[200100], bit2[200100];
void cng(int x, modint w) {
    for (; x<=n; x+=x&-x) {
        bit[x] += w;
    }
}
modint qry(int x) {
    modint res = 0;
    for (; x; x-=x&-x) {
        res += bit[x];
    }
    return res;
}
void cng2(int x, modint w) {
    for (; x<=n; x+=x&-x) {
        bit2[x] += w;
    }
}
modint qry2(int x) {
    modint res = 0;
    for (; x; x-=x&-x) {
        res += bit2[x];
    }
    return res;
}
void solve() {
    cin >> n;
    bastr<int> dsc;
    for (int i=1; i<=n; ++i) {
        cin >> a[i];
        dsc += a[i];
    }
    sort(dsc.begin(), dsc.end());
    dsc.erase(unique(dsc.begin(), dsc.end()), dsc.end());
    for (int i=1; i<=n; ++i) {
        a[i] = lower_bound(dsc.begin(), dsc.end(), a[i])-dsc.begin()+1;
    }
    iota(id+1, id+n+1, 1);
    sort(id+1, id+n+1, [&](int x, int y) {
        return a[x] < a[y];
    });
    for (int i=1,now=n; i<=n; ++i) {
        while (now >= 1 && a[now] <= a[id[i]]) {
            now--;
        }
        pos[id[i]] = now;
    }
    for (int i=1; i<=n; ++i) {
        bit[i] = 0;
    }
    for (int i=1; i<=n; ++i) {
        dp[i] = qry(a[i]-1)+1;
        cng(a[i], dp[i]);
    }
    for (int i=1; i<=n; ++i) {
        bit[i] = bit2[i] = 0;
    }
    modint ans = 0;
    for (int i=n; i>=1; --i) {
        dp2[i] = qry(cs(dsc)-a[i])+1;
        cng(cs(dsc)-a[i]+1, dp2[i]);
        if (pos[i] > i) {
            dp3[i] = qry(cs(dsc)-a[pos[i]]+1)+qry2(a[pos[i]]-1)-qry2(a[i]);
            cng2(a[i], dp3[i]);
        }
        else {
            dp3[i] = 0;
        }
        ans += dp[i]*dp3[i];
    }
    cout << ans.val() << "\n";
}

link。

容易发现,如果将 $x$ 写作 $\displaystyle \prod_{i = 1}^k p_i^{\alpha_i}$ 的形式,$\displaystyle J(x) = 1+\sum p_i^{\alpha_i}+\sum\sum p_i^{\alpha_i}p_j^{\alpha_j}+\dots = \sum_{T \in 2^S} \sum_{i \in T} p_i^{\alpha_i}$,其中 $S = \{1, \dots, k\}$。写到这里可以发现 Joker function 是个 multiplicative function,所以 Joker function 又可以写作 $\displaystyle J(x) = \prod_{i = 1}^k J(p_i^{\alpha_i}) = \prod_{i = 1}^k \left(p_i^{\alpha_i}+1\right)$。

考虑 dp,设 $f[i][j]$ 表示用前 $i$ 个因数通过上述形式凑出 $j$ 的方案数,转移很平凡,具体看代码。这个 dp 可以用数据结构来存。这个 dp 复杂度的保证在于 $\displaystyle \max_{i \in [1, 10^{12}]} {\sigma_0(i)} = 6720$,这就叫人比较无语。至于判断一个数是否是 $p^k+1$ 的形式,可以根号分治来做,记阈值为 $T = 10^6$,对于 $\leqslant T$ 的元素,可以线性筛出所有质数,标记所有质数的次幂,数量级粗略是 $O(T\log T)$,完全跑不满;对于 $> T$ 的元素,因为 $(T+1)^2 > 10^{12}$,所以只需要判断是否为质数即可,Miller Rabin 或者其他 nb 的判素数方法即可。写不来 MR 所以直接蒯了个 mrsrz 的。

namespace Miller_Rabin{
    const int P[]={2,3,7,97,19260817};
    inline LL mul(LL a,LL b,const LL&md){
        LL tmp=a*b-(LL)((long double)a*b/md+.5)*md;
        return tmp+(tmp>>63&md);
    } 
    inline LL pow(LL a,LL b,const LL&md){
        LL ret=1;
        for(;b;b>>=1,a=mul(a,a,md))if(b&1)ret=mul(ret,a,md);
        return ret;
    }
    bool test(LL a,LL p){
        LL x=p-1;
        int d=0;
        while(!(x&1))++d,x>>=1;
        LL t=pow(a,x,p);
        while(d--){
            const LL ls=t;
            t=mul(t,t,p);
            if(t==1&&ls!=1&&ls!=p-1)return 0;
        }
        return t==1;
    }
    bool check(LL n){
        if(n<2)return 0;
        if(n==2||n==3||n==7||n==97||n==P[4])return 1;
        for(int i=0;i<5;++i)if(n%P[i]==0)return 0;
        for(int i=0;i<5;++i)
        if(!test(P[i]%n,n))return 0;
        return 1;
    }
}
using Miller_Rabin::check;
LL A;
bitset<1000100> tag;
int prm[1000100], cnt;
map<LL, int> mp; // a in it iff a = p^k
bastr<LL> offset[2000100];
void sieve(int n) {
    for (int i=2; i<=n; ++i) {
        if (!tag.test(i)) {
            prm[++cnt] = i;
        }
        for (int j=1; j<=cnt && i<=n/prm[j]; ++j) {
            tag.set(i*prm[j]);
            if (i%prm[j] == 0) {
                break;
            }
        }
    }
    for (int i=1; i<=cnt; ++i) {
        for (LL j=prm[i]; j<=1e12; j*=prm[i]) {
            mp[j] = i;
        }
    }
}
void executer() {
    cin >> A;
    map<LL, LL> dp[2];
    bastr<int> eff;
    for (LL i=1; i<=A/i; ++i) {
        if (A%i == 0) {
            if (mp.count(i-1)) {
                offset[mp[i-1]] += i;
                eff += mp[i-1];
            }
            if (i*i != A) {
                if (mp.count(A/i-1)) {
                    offset[mp[A/i-1]] += A/i;
                    eff += mp[A/i-1];
                }
                else if (check(A/i-1)) {
                    offset[++cnt] += A/i;
                    eff += cnt;
                }
            }
        }
    }
    sort(eff.begin(), eff.end());
    eff.erase(unique(eff.begin(), eff.end()), eff.end());
    dp[0][1] = 1;
    int cur = 1;
    for (auto u : eff) {
        dp[cur] = dp[cur^1];
        for (auto v : dp[cur^1]) {
            for (auto p : offset[u]) {
                if (A/p >= v.first && A%(v.first*p) == 0) {
                    dp[cur][v.first*p] += v.second;
                }
            }
        }
        cur ^= 1;
    }
    cout << dp[cs(eff)&1][A] << "\n";
}

link。

还算萌,但是代码有些难写……

你首先会想要

int n, m, fa[20][300100], pa[300100], dep[300100], cnt[900100];
int ldf[300100], rdf[300100], dfc, qwq;
vi<int> G[300100];
struct path {
    int x, y, sx, sy, lca;
};
void dfs(int x, int pre) {
    fa[0][x] = pa[x] = pre;
    dep[x] = dep[pre]+1;
    ldf[x] = ++dfc;
    for (int i=1; i<20; ++i) {
        fa[i][x] = fa[i-1][fa[i-1][x]];
    }
    for (auto y : G[x]) {
        if (y != pre) {
            dfs(y, x);
        }
    }
    rdf[x] = dfc;
}
int lca(int x, int y) {
    if (dep[x] < dep[y]) {
        swap(x, y);
    }
    for (int i=19; i>=0; --i) {
        if (dep[fa[i][x]] >= dep[y]) {
            x = fa[i][x];
        }
    }
    if (x == y) {
        return x;
    }
    for (int i=19; i>=0; --i) {
        if (fa[i][x] != fa[i][y]) {
            x = fa[i][x], y = fa[i][y];
        }
    }
    return pa[x];
}
int jump(int x, int d) {
    if (d < 0) {
        return n+(++qwq);
    }
    for (int i=19; i>=0; --i) {
        if ((d>>i)&1) {
            x = fa[i][x];
        }
    }
    return x;
}
int bit[300100];
void add(int x, int y) {
    for (; x<=n; x+=x&-x) {
        bit[x] += y;
    }
}
int qry(int x) {
    int res = 0;
    for (; x; x-=x&-x) {
        res += bit[x];
    }
    return res;
}
int qry(int l, int r) {
    return qry(r)-qry(l-1);
}
void executer() {
    cin >> n;
    for (int i=1,x,y; i<n; ++i) {
        cin >> x >> y;
        G[x].push_back(y);
        G[y].push_back(x);
    }
    dfs(1, 0);
    cin >> m;
    vi<path> paths(m);
    for (auto& it : paths) {
        cin >> it.x >> it.y;
        it.lca = lca(it.x, it.y);
        it.sx = jump(it.x, dep[it.x]-dep[it.lca]-1);
        it.sy = jump(it.y, dep[it.y]-dep[it.lca]-1);
        if (it.sx > it.sy) {
            swap(it.sx, it.sy), swap(it.x, it.y);
        }
    }
    sort(paths.begin(), paths.end(), [&](const path& lhs, const path& rhs) {
        return dep[lhs.lca] < dep[rhs.lca]
               || (dep[lhs.lca] == dep[rhs.lca] && lhs.lca < rhs.lca)
               || (lhs.lca == rhs.lca && lhs.sx > rhs.sx)
               || (lhs.sx == rhs.sx && lhs.sy > rhs.sy);
    });
    LL ans = 0;
    for (int i=0; i<m;) {
        int j = i;
        while (j < m-1 && paths[j+1].lca == paths[j].lca) {
            j++;
        }
        LL now = 0;
        for (int l=i; l<=j;) {
            int r = l;
            while (r < j && paths[r+1].sx == paths[r].sx) {
                r++;
            }
            for (int k=l; k<=r; ++k) {
                ans += now-cnt[paths[k].sy];
            }
            for (int k=l; k<=r; ++k) {
                cnt[paths[k].sx]++, cnt[paths[k].sy]++;
            }
            now += r-l+1, l = r+1;
        }
        for (int k=i; k<=j; ++k) {
            cnt[paths[k].sx] = cnt[paths[k].sy] = 0;
        }
        for (int k=i; k<=j; ++k) {
            ans += qry(ldf[paths[k].lca], rdf[paths[k].lca]);
            if (paths[k].sx <= n) {
                ans -= qry(ldf[paths[k].sx], rdf[paths[k].sx]);
            }
            if (paths[k].sy <= n) {
                ans -= qry(ldf[paths[k].sy], rdf[paths[k].sy]);
            }
        }
        for (int k=i; k<=j; ++k) {
            add(ldf[paths[k].x], 1), add(ldf[paths[k].y], 1);
        }
        i = j+1;
    }
    cout << ans << "\n";
}