莫队

發表於 2020-08-28 更新於 2025-06-05 分類於 algorithm Disqus：文章字數： 714 所需閱讀時間 ≈ 2 分鐘

题目

啊被牛客题目描述坑了，把关键信息藏在数据范围里，不过顺便学了一下莫队。

牛客题目https://ac.nowcoder.com/acm/contest/7079/A ，知道真实题意以后是个水题

大意

指莫队要解决的问题不是牛客的题。

给长n的数组，没有改动，q次询问区间的不同值的个数

q和n都1e5

解法

我看到mex都很慌，这玩意虽然看上去是个质朴的东西，但涉及到的都不太好搞

上面直接暴力，就 O(q n)。

莫队实际上就是优雅的分块，然后再暴力。

离线请求。
根据请求的[左端点/sqrt(n), 右端点]排序
暴力求解

意味着，每次对于 左端点/sqrt(n)相同的时候，时间复杂度= sqrt(n)*询问个数 + n

总时间复杂度 sqrt(n)*q + n*sqrt(n)

一些优化，通过奇偶来决定右端点是正序还是逆序列。

代码

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
#define ten5 100000+10
#define rep(i,a,n) for (int i=a;i<n;i++)
#define iif(c,t,f) ((c)?(t):(f))
#define per(i,a,n) for (int i=n-1;i>=a;i--)


const int MAXN = 30005; // 数组大小
const int MAXQ = 200005; // 询问次数
const int MAXM = 1000005; // 值记录

int sq; // sqrt(n);
struct query{ // 把询问以结构体方式保存
  int l, r, id;
  // 不是按照 l 第一序列，r第二序
  // 而是 l/sq 第一序，r第二序列，对于相同的l/sq 其中根据l/sq的奇偶性来决定r是正向序还是逆向序
  bool operator<(const query &o) const { // 重载<运算符，奇偶化排序
    // 这里只需要知道每个元素归属哪个块，而块的大小都是sqrt(n)，所以可以直接用l/sq
    if (l / sq != o.l / sq)
      return l < o.l;
    if (l / sq & 1)
      return r < o.r;
    return r > o.r;
  }
} Q[MAXQ];
int A[MAXN];// 原始数组
int ans[MAXQ]; // 离线结果
int Cnt[MAXM];
int cur;// 总计
int l = 1, r = 0; // 左右游标

void add(int p) {
  if (Cnt[A[p]] == 0)
    cur++;
  Cnt[A[p]]++;
}
void del(int p) {
  Cnt[A[p]]--;
  if (Cnt[A[p]] == 0)
    cur--;
}
ll read(){
  ll r;
  scanf("%lld",&r);
  return r;
}
int main() {
  int n = read();
  sq = sqrt(n);
  rep(i,1,n+1){
    A[i] = read();
  }
  int q = read();
  rep(i,0,q){
    Q[i].l = read();
    Q[i].r = read();
    Q[i].id = i; // 把询问离线下来
  }
  sort(Q, Q + q); // 排序
  // 每个 按 l/sq的区块分的区间内，r的总变化最多n, l的变化 = sqrt(n) * 个数
  // 总代价 = sum 所有区块的代价 = 区块数 * r的每个区块变化 + sqrt(n) * 总个数 = n^(3/2) + q * n^(1/2);
  rep(i,0,q){
    while (l < Q[i].l)
      del(l++);
    while (l > Q[i].l)
      add(--l);
    while (r < Q[i].r)
      add(++r);
    while (r > Q[i].r)
      del(r--);
    ans[Q[i].id] = cur; // 储存答案
  }
  rep(i,0,q){
    printf("%d\n", ans[i]); // 按编号顺序输出
  }
  return 0;
}

参考

各种莫队博文

project euler 101 Optimum polynomial

發表於 2020-07-01 更新於 2025-06-05 分類於 Project Euler Disqus：文章字數： 953 所需閱讀時間 ≈ 3 分鐘

题目

https://projecteuler.net/problem=101

大意

给你一个多项式表达式

然后只对前面几项进行最小幂次拟合，拟合的函数第一个和原函数不同的输出值的和

原题目很简单，11x11的Vandermonde矩阵求个逆矩阵，甚至都不用处理分数就搞完了

hackerrank

然后我发现了这个

https://www.hackerrank.com/contests/projecteuler/challenges/euler101/problem

直接把矩阵拉到了 3000x3000,然后多项式的系数是传入的

首先Vandermonde就需要一个高效求逆的办法，最后再O(n^2)算算

高效求逆?

为什么要, 因为我们有了原表达式，要去做拟合实际就是

${\begin{pmatrix}
{a_0}\\
{a_1}\\
{a_2}\\
{\vdots}\\
{a_{n-1}}\\
\end{pmatrix}
}^T\begin{pmatrix}
{1}&{1}&{\cdots}&{1}\\
{x_{0}}&{x_{1}}&{\cdots}&{x_{n-1}}\\
{x_{0}^2}&{x_{1}^2}&{\cdots}&{x_{n-1}^2}\\
{\vdots}&{\vdots}&{}&{\vdots}\\
{x_{0}^{n-1}}&{x_{1}^{n-1}}&{\cdots}&{x_{n-1}^{n-1}}\\
\end{pmatrix}={
\begin{pmatrix}
{y_0}\\
{y_1}\\
{y_2}\\
{\vdots}\\
{y_{n-1}}
\end{pmatrix}}^T$

要去计算左侧的a,如果能快速求逆，那么有

${\begin{pmatrix}
{a_0}\\
{a_1}\\
{a_2}\\
{\vdots}\\
{a_{n-1}}\\
\end{pmatrix}}^T
={\begin{pmatrix}
{y_0}\\
{y_1}\\
{y_2}\\
{\vdots}\\
{y_{n-1}}
\end{pmatrix}}^T
\begin{pmatrix}
{1}&{1}&{\cdots}&{1}\\
{x_{0}}&{x_{1}}&{\cdots}&{x_{n-1}}\\
{x_{0}^2}&{x_{1}^2}&{\cdots}&{x_{n-1}^2}\\
{\vdots}&{\vdots}&{}&{\vdots}\\
{x_{0}^{n-1}}&{x_{1}^{n-1}}&{\cdots}&{x_{n-1}^{n-1}}\\
\end{pmatrix}^{-1}$

也就是求完逆以后，只需要O(n^2)来得到a

$V(x_0,x_1,\cdots ,x_{n-1})=\begin{bmatrix} {1}&{1}&{\cdots}&{1} \\
{x_{0}}&{x_{1}}&{\cdots}&{x_{n-1}}\\
{x_{0}^2}&{x_{1}^2}&{\cdots}&{x_{n-1}^2}\\
{\vdots}&{\vdots}&{}&{\vdots}\\
{x_{0}^{n-1}}&{x_{1}^{n-1}}&{\cdots}&{x_{n-1}^{n-1}}\\ \end{bmatrix}$

有 $ V(x_0,x_1,\cdots ,x_{n-1})=\prod _{n > i > j \geq 0}(x _{i}-x _{j})$

众所周知 $A^{-1} = \frac{1}{|A|} A^*$,然而这暴力算是4次方复杂度，似乎比玩高斯消元3次方还要久

换个方法拉格朗日插值

想法很简单比如过点 (1,1) (4,7) (9,100)的二次函数

直接表达式$f(x) = k_1(x-4)(x-9)+k_2(x-1)(x-9) + k_3(x-1)(x-4)$

注意到x取其中一个点时，求和的表达式只有一个不为0

也很明显 $k_i = y_i\prod\limits_{j\not =i}\dfrac{1}{x_i-x_j}$

$f(x)=\sum\limits_{i=0}^ky_i\prod\limits_{j\not =i}\dfrac{x-x_j}{x_i-x_j}$, 也可以知道它的逆矩阵就是求和每一部分的展开式的系数，(令 $y = [0 … 0 1 0 … 0]$ 即可知

回到题目$x_i = i+1, y_i = \sum\limits_{p=0}^n A_p(i+1)^p$

$f(x)=\sum\limits_{i=0}^k[y_i\prod\limits_{j\not =i}\dfrac{x-j-1}{i-j}]$

注意到我们其实只需要求 $f(k+1)$ , 计算所有$y_i$ O(n^2), 但是右侧看似复杂的$x_i$乘积，因为这里的取值，变成了只需要做阶乘运算即可。可类似前缀和。再配合一点乘法逆元即可

代码

N = 0
A = []
inv = []
po = []
poinv = []
MOD = 1000000007


def u(v):
    r = 0
    vi = 1
    for i in range(0, N+1):
        r += A[i] * vi
        r %= MOD
        vi *= v
        vi %= MOD
    return r


def init():
    # init inv
    global inv
    inv = []
    inv.append(0)
    inv.append(1)
    for i in range(2, N+10):
        inv.append(((MOD - MOD // i) * inv[MOD%i]) % MOD)
    # init power
    global po
    po = []
    po.append(1)
    for i in range(1, N+10):
        po.append( (po[i-1] * i) % MOD )

    global poinv
    poinv = []
    poinv.append(1)
    for i in range(1, N+10):
        poinv.append( (poinv[i-1] * inv[i]) % MOD )


def fitn(y, n):
    r = 0
    for i in range(1, n+1):
        mul = y[i-1]
        # 分子
        mul *= po[n] * inv[n+1-i]
        # 分母乘法逆元
        mul *= poinv[i-1] * poinv[n-i]
        if (n - i) % 2 == 1:
            mul *= -1
        mul %= MOD
        r += mul
        r %= MOD
    return (r+MOD)%MOD


def work():
    init()
    ans = 0
    y = []
    for v in range(1, N+1):
        arr.append(u(v))

    for i in range(1, N+1):
        r = fitn(y, i)
        print(r, end=' ')
        ans += r
    return ans


def pe():
    global N, A
    N = 10
    A = [1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1]
    print(work())


def main():
    global N, A
    N = int(input())
    A = list(map(int, input().split()))
    work()


# pe()
main()

至此我们通过了较难版本的hackerrank 上的pe101

参考

http://www.vesnik.math.rs/vol/mv19303.pdf

https://proofwiki.org/wiki/Inverse_of_Vandermonde_Matrix

project euler 84 Monopoly odds 不会了算概率是认真的吗？

發表於 2020-06-27 更新於 2025-06-05 分類於 Project Euler Disqus：文章字數： 332 所需閱讀時間 ≈ 1 分鐘

题目

https://projecteuler.net/problem=84

大意

玩大富翁

有一些位置会触发跳转（跳转逻辑贼复杂）

问落在哪三个位置的概率最大

题解

我都不想贴代码，这玩意真是靠大量数据“算”概率的咯。

所以伪随机数还是能作为概率意义上的“随机数”是因为题目假设赛子均匀，而随机数分布的期望是均匀，数据量大了就符合统计逻辑？

想说有不少是模拟行走去记录次数然后处以总步数，然而这样的操作如果把4面换到6面很大才能收敛到题目已经告诉你的102400

所以虽然4面可以比较小的数就得到答案，但感觉更多是因为小，所以“猜”对的。

感觉还是应该所有点 2.5%，然后每一轮进行所有格子走一次，直接计算概率转移，而不是模拟random，从而直到收敛（这也是能判断的了，不像上面计次难以判断是否收敛，只能靠设置大，更大，再大的步数）。

疑惑的是，很多题目都是给小数据，要求大数据，怎么这一道题是给6求4。

其实读题,题太长才是我觉得的问题。XD

project euler 100 Arranged probability

發表於 2020-06-23 更新於 2025-06-05 分類於 Project Euler Disqus：文章字數： 2k 所需閱讀時間 ≈ 7 分鐘

题目

https://projecteuler.net/problem=100

大意

$\frac{x(x-1)}{y(y-1)} = \frac{1}{2}$,x,y均为正整数，求首个满足表达式且$x>10^{12}$的x的值

读懂本文你可能需要的前置知识

pell方程和连分数 project euler 66

解

显然通过等式变换的基本操作有

$(2y-1)^2-2(2x-1)^2 = -1$

这不就是 project euler的66题很像的pell 方程吗！

一点比较好的是$10^{12}$,其实可以双指针暴搜,O(x)时间复杂度，O(1)空间暴力算几个小时算出来,但当然如果是暴力也就不会写这篇文章了

稍加搜索也能搜到这个 https://www.alpertron.com.ar/QUAD.HTM 也就是说$ax^2+bxy+cy^2+dx+ey+f=0$的等式都能解，

当然直接用网页也不需要这篇文章了，本文还是以证明为核心

简化

考虑 $ax^2+bxy+cy^2+dx+ey+f=0$ 二元二次方程(以下过程中所有字母均为整数)

因为对称性，我们只考虑$a \neq 0$，否则交换x和y。

不然就是形如$axy+bx+cy+d=0$的方程了，而这类方程总能化为$(ax+b)(cy+d) = e$的形式，变为拆分e的。众所周知比较容易，本文就不涉及细节了。

$4a(ax^2+bxy+cy^2+dx+ey+f)=0$

$4a^2x^2+4abxy+4acy^2+4adx+4aey+4af=0$

$(2ax+by+d)^2 +(4ac-b^2)y^2+(4ae-2bd)y+(4af-d^2)= 0$

这形如

$x^2+ay^2+by+c=0$, 如果这里$y^2$的系数$a=0$,那么也是$x^2+ay+b=0$的简单情况，讨论$x \in [0,a)$，这里不细聊了,以下只考虑$a\neq 0$

$4a^2(x^2+ay^2+by+c) = 0$

$(2ax)^2+4a^3y^2+4a^2by+4a^2c = 0$

$(2ax)^2+a(2ay+b)^2+(4a^2c-ab^2) = 0$

这形如

$x^2+ay^2+c=0$

最终又回到了这个式子，如果这里的a是正数，显然是一个椭圆方程，有限范围内的一些尝试,不是本次讨论的，暂时不展开了

那么最终回到了

$x^2-dy^2 = k$,其中$d > 0,k \neq 0$, $d$为非平方数,否则是平方数间隔问题

温故知新

看看之前project euler 66 的证明过程

pell方程基础解
pell方程解存在性
连分数分子分母递推式
连分数，分子分母相关性质(分子分母互质性,原无理数与渐进数差值表达式，差值绝对值单调递减性，连分数是最给定分母以及更小分母中最接近原无理数，任何满足“某不等式”的数字是渐进数（这个不等式是pell方程的必要条件）)
Lagrange’s Theorem 二次无理数与有循环节的连分数
Galois’ Theorem 纯循环连分数与reduced二次有理数
根号n的连分数性质（精确的$\frac{P_n+\sqrt{d}}{Q_n}$表示和形如pell方程与$Q_n$的等式）
分类讨论$Q_n$与n的取值最后的结论

对比有

1&2不再是一定存在的，我们会讨论如果存在，那么具有性质和基础解

3只涉及连分数的性质，之前证明可复用

4同样只涉及连分数和不等式，最终结论任何最简分数$\frac{h}{k}$如果满足 $|\sqrt{d} - \frac{h}{k}| < \frac{1}{2k^2}$ 那么它是个渐进分数,同样复用

5,6都是连分数循环性质，reduced性质和二次无理数的关系，可复用

7是根号n的连分数表示性质，和$h_{n-1}^2-dk_{n-1}^2=(-1)^nQ_n$的证明，可复用

8再回到$x^2-dy^2=k$

如果存在

若存在$(x_0,y_0),x_0 > 0,y_0 > 0$使得 $x_0^2-dy_0^2 = k$, 因为有最小的$(x_1,y_1)$满足$x_1^2-dy_1^2=1$ (证明见pe 66的证明)

那么通过$(x_0+y_0\sqrt{d})(x_1+y_1\sqrt{d})$能得到新的解

$(x_0^2-dy_0^2)(x_1^2-dy_1^2) = (x_0x_1+dy_0y_1)^2 - d(x_0y_1+x_1y_0)^2 = k$，所以如果有解那么有无限组解

设$(x_i,y_i)$是满足$x^2-dy^2 = k$的最小解,$(x_j,y_j)$是另一个满足方程的解,有$x_i<x_j,y_i<y_j$

若$(x_j,y_j)$不是由最小的$(x_i,y_i)$和$(x_1,y_1)$生成的,则

$(x_i+y_i\sqrt{d})(x_1+y_1\sqrt{d})^n < (x_j+y_j\sqrt{d}) < (x_i+y_i\sqrt{d})(x_1+y_1\sqrt{d})^{n+1}$

同时乘$(x_1-y_1\sqrt{d})^n$

$x_i+y_i\sqrt{d} < (x_j+y_j\sqrt{d})(x_1-y_1\sqrt{d})^n < (x_i+y_i\sqrt{d})(x_1+y_1\sqrt{d})$

因为中间的也满足 $x+y\sqrt{d}$的形式

若$(x_j+y_j\sqrt{d})(x_1-y_1\sqrt{d})^n$ 也是一组生成的解

令$(x_j+y_j\sqrt{d})(x_1-y_1\sqrt{d})^n = x_p+y_p\sqrt{d}$

有$(x_j-y_j\sqrt{d})(x_1+y_1\sqrt{d})^n = x_p-y_p\sqrt{d}$

** TODO 不会证明了失败 **

证明

$1=\frac{x_0^2-dy_0^2}{x_1^2-dy_1^2} $

$= \frac{x_0+y_0\sqrt{d}}{x_1+y_1\sqrt{d}} \cdot \frac{x_0-y_0\sqrt{d}}{x_1-y_1\sqrt{d}}$

$= \frac{(x_0+y_0\sqrt{d})(x_1-y_1\sqrt{d})}{-1} \cdot \frac{(x_0-y_0\sqrt{d})(x_1+y_1\sqrt{d})}{-1}$

$= ((x_0x_1-y_0y_1d)+(x_1y_0-x_0y_1)\sqrt{d})((x_0x_1-y_0y_1d)-(x_1y_0-x_0y_1)\sqrt{d})$

$= (x_0x_1-y_0y_1d)^2-(x_1y_0-x_0y_1)^2 d $

说明两个不同的-1的解一定是通过乘一个=1的解得到的

说明存在一个=-1的基础解，其它解是靠基础解乘上(=1的基础解的幂次)得到的

主要关心 = -1

同时我们在66证明过 $h_{n-1}^2-dk_{n-1}^2 = (-1)^nQ_n$

**注记：最开始66那一版本,的渐进分数的部分直接使用的是pell方程中$\sqrt{d}$,其中d是非平方正整数，但下面的部分证明和证明实际上只需要保证其中部分的参数只要是无理数即可。所以进行了修改 **

若$x^2-\alpha^2 y^2 = \beta$, 且存在解，若$|\beta| < \alpha $,则$\frac{x}{y}$ 一定是 $\alpha$的渐近分数。

证明

若$\beta > 0$,有 $x^2-\alpha^2y^2 > 0 $即$x>y\alpha$

有$|\frac{x}{y} - \alpha| = \frac{\beta}{y(x+y\alpha)} < \frac{\beta}{2y^2\alpha} < \frac{1}{2y^2}$

前面证明过(66的6.2.3.8)满足这个表达式的对应的最简分数$\frac{x}{y}$ 一定是渐近连分数

若$\beta < 0$

$y^2 - \frac{1}{\alpha^2}x^2 = \frac{-\beta}{\alpha^2}$,注意到$|\frac{-\beta}{\alpha^2}| < \frac{1}{ \alpha } $ 也就回到上面的情况

综上，证毕。

以上的证明说明了如果$x^2-dy^2=-1$如果有解,那么一定是一个最简分数

同样根据66的过程我们连分类讨论的内容都能复用

$h_{n-1}^2-dk_{n-1}^2=(-1)^nQ_n = (-1)^n, n = kl $,l为$\sqrt{d}$的最小周期长度,

也就是再看个奇偶，最后得到了pell方程和连分数比较紧密的结论

$l$长度为偶数无解

$l$长度为奇数, $(x,y) = {h_{(2m+1)l-1},k_{(2m+1)l-1}}$,m为非负整数

回到题目

$x^2-2y^2 = -1$

关于$\sqrt{2}$的连分数展开众所周知了也(66连分数开篇就是)

从前面得到的结论看也可以通过$(x_0+\sqrt{d}y_0)(x_1+\sqrt{d}y_1)^n$来生成

$(1+\sqrt{2})(3+2\sqrt{2})^n$

$(x+y\sqrt{2})(3+2\sqrt{2})^n$

$(x_n,y_n) = (3x_{n-1}+4y_{n-1},2x_{n-1}+3y_{n-1})$

代码

嗯高亮有点问题整除识别成注释了

x = 1
y = 1

while True:
    x,y = 3*x+4*y,2*x+3*y;
    print(x,y);
    if (x+1)//2 > 1000000000000:
        print((y+1)//2)
        break

参考

https://oeis.org/A031396

https://en.wikipedia.org/wiki/Pell%27s_equation#The_negative_Pell_equation

https://mathworld.wolfram.com/PellEquation.html

总结&感受

有搜到有人通过枚举前面的小解，发现了增长比例有趋近于定值的趋势。哎我连找规律也不会了吗

没聊到其它非1和非-1时的情况

似乎$x^2-34y^2=-1$不存在解,所以可能无法证明它一定有解，但是如果有解，上面证明了一定是通过连分数的渐进分数找到的,也可以通过基础解生成

project euler 078 Coin partitions

發表於 2020-06-10 更新於 2025-06-05 分類於 Project Euler Disqus：文章字數： 1.9k 所需閱讀時間 ≈ 6 分鐘

题目

https://projecteuler.net/problem=78

4 = 4 = 3+1 = 2+2 = 2+1+1=1+1+1+1,p(4)=5

5 = 5 = 4+1 = 3+2 = 3+1+1 = 2+2+1 = 2+1+1+1 = 1+1+1+1+1,p(5)=7

求最小的n使得 p(n) % 1_000_000 = 0

rust

dp递推暴力似乎并不可解,警告下面的代码只要没找到答案就会一直增加vector大小，所以建议尽早掐断，不然会吃很大内存XD,可能造成系统卡死或者假卡死

use std::cmp;


// dp[n][x] + (n-x <= x?) + sum(dp[n-x][1 -> min(x,n-x-1)])


fn main() {
    let mut dp = Vec::new();
    let mut dpsum = Vec::new();
    dp.push(Vec::new());
    dpsum.push(Vec::new());
    let mut i = 1;
    loop{
        dp.push(Vec::new());
        dpsum.push(Vec::new());
        dp[i].push(0);
        dpsum[i].push(0);
        let mut cnt = 0;
        for j in 1..i+1 {
            dp[i].push(0);
            dpsum[i].push(0);
            if i == j {
                dp[i][j] = 1;
            }else{
                dp[i][j] = dpsum[i-j][cmp::min(j,i-j)];
            }
            cnt += dp[i][j];
            cnt %= 1_000_000;
            dpsum[i][j] = (dpsum[i][j-1]+dp[i][j]) % 1_000_000;
        }
        println!("{}=>{}",i,cnt);
        if cnt % 1_000_000 == 0 {
            println!("{}=>{}",i,cnt);
            return ;
        }
        i+=1;
    }
}

题解

分拆

实际上可以把n看成，拆了$a_1$个1,$a_2$个2…

所以$p(n)$可以看成 $n = 1a_1+2a_2+3a_3+…+na_n$的解 $(a_1,a_2,…,a_n)$的个数

也是 $(1+x+x^2+x^3+…)(1+x^2+x^4+x^6+…)(1+x^3+x^6+x^9+…)…$ 结果中 $x^n$的系数,因为相当于从每个括号中选取了一项假设第i个括号选取的是$x^{ij}$ 那么对应上面 $a_i = j$，也对应到相关的解

也就是 $\sum_{n=0}^{\infty} p(n)x^n = (1+x+x^2+…)(1+x^2+x^4+…)… = \frac{1}{1-x} \frac{1}{1-x^2} … = \prod_{n=1}^{\infty} \frac{1}{1-x^n} $

互异分拆

如果分拆结果要求不能出现相同的数字(互异分拆 partitions into distinct parts)，则相当于上面$a_i$增加限制$a_i\in\{0,1\}$

也就是 $\sum_{n=0}^{\infty} q(n)x^n = (1+x)(1+x^2)… = \prod_{n=1}^{\infty} (1+x^n) $

$p(n) = $拆分成偶数个互异的数的方案数 + 拆分成奇数个互异的数方案数

考虑一个互异分拆$n=n_1+n_2+…+n_s$ (注意这里是$n_i$不是$a_i$,实际是$a_{n_i} = 1$的那些)

$x^{n} = x^{n_1}x^{n_2}…x^{n_s}$

$(-1)^sx^{n} = (-x^{n_1})(-x^{n_2})…(-x^{n_s})$

考虑 $\sum_{n=0}^{\infty} g(n)x^n = (1-x)(1-x^2)… = \prod_{n=1}^{\infty} (1-x^n) $ 从值的构成上它依然表示的是互异分拆

注意到:$a_{n_i} = 0$在乘积贡献上是1,$a_{n_i}=1$在乘积贡献上是$-1$,所以上述每个互异分拆对$g(n)$的贡献是 $(-1)^s$

也就是$g(n) = $拆分成偶数个互异的数方案数 - 拆分成奇数个互异的数方案数

g(n)

下面我们要看看欧拉大佬的神之推导

先上 Ferrer 图

比如 $ 11 = 5+4+2$

1
2
3

xxxxx
xxxx
xx

只考虑互异分拆，把底部的长度记为b(上图是1),从最大递减1的长度记为s(上图是5,4 所以s=2)

如果一个图 $s \le b$,那么可以把最下面一排移到右侧,上图变为

1 2	xxxxx x xxxx x

注意性质

移动以后，因为少了一行，分拆的奇偶性变了。
b至少增加1, s = 原来的b ，所以新的图是 s>b的。

仅有一种情况是不能移动的（因为要移动的最后一行点数比移动后剩余的行数大1）如下,注意到需要s排满和s排满，所以对于分拆的一个n，至多出现一次，学过等差数列求和的也可以看出满足这种的n的公式。

$n = \frac{s(s+(2s-1))}{2}$

1
2
3

xxxxx
xxxx
xxx

以上说明除了这特殊的以外，所有$s \le b$的都有唯一对应的 $s > b$的互异ferrer图

再看 $s>b$的,例如(s=3,b=2)

1
2
3

xxxxxx
xxxxx
xxx

我们把右侧45度的点移到下面一排，有

xxxxx
xxxx
xxx
xx

同样有性质

移动以后因为多了一行，分拆的奇偶性变了
新的s >=旧的s，新的b = 旧的s，所以有新的s>=新的b

同样有特殊的$s+1=b$的排满时，不能移动（否则移动后最下两行点数相等）

1
2
3

xxxxxx
xxxxx
xxxx

同样学过等差数列求和的，对于一个n最多有一个，n的表达式也呼之欲出了

$n = \frac{s((s+1)+2s)}{2}$

综上除了两种特殊的情况，其它的都能找到一一对应的，所以偶奇相减除了特殊的外结果为0

所以就有神奇的展开

$\prod_{n=1}^{\infty}(1-x^n) = 1+\sum_{s=1}^{\infty} (-1)^s(x^{\frac{3s^2-s}{2}}+x^{\frac{3s^2+s}{2}})$

回到目标p(n)

$(\sum_{n=0}^{\infty} p(n)x^n)(\cdot \prod_{n=1}^{\infty} 1-x^n) = (\prod_{n=1}^{\infty} )\cdot(\frac{1}{1-x^n} \cdot \prod_{n=1}^{\infty} 1-x^n) = 1 $ (…真的能这样乘吗。。

$(\sum_{n=0}^{\infty} p(n)x^n)(1+\sum_{s=1}^{\infty} (-1)^s(x^{\frac{3s^2-s}{2}}+x^{\frac{3s^2+s}{2}})) = 1 $

也就意味着左边乘以后对于$n>0$ 乘积后的系数都是$a_n = 0$

$a_n x^n = (p(n)x^n\cdot 1 + \sum_{s=1}^{\frac{3s^2-s}{2} \le n}[ p(n-\frac{3s^2-s}{2}) x^{n-\frac{3s^2-s}{2}} (-1)^s \cdot x^{\frac{3s^2-s}{2}} ]+ \sum_{s=1}^{\frac{3s^2+s}{2} \le n}[ p(n-\frac{3s^2+s}{2}) x^{n-\frac{3s^2+s}{2}} (-1)^s \cdot x^{\frac{3s^2+s}{2}} ])$

$0=a_n = p(n) + \sum_{s=1}^{\frac{3s^2-s}{2} \le n}[ p(n-\frac{3s^2-s}{2}) (-1)^s ]+ \sum_{s=1}^{\frac{3s^2+s}{2} \le n}[ p(n-\frac{3s^2+s}{2}) (-1)^s ])$

由此我们有了$p(n)$的递推式,时间复杂度像是$n^1.5$(因为 s与n的不等式是s平方增长，空间就只需要$O(n)$

$p(n) = \sum_{s=1}^{\frac{3s^2-s}{2} \le n}[(-1)^{s+1} p(n-\frac{3s^2-s}{2})]+ \sum_{s=1}^{\frac{3s^2+s}{2} \le n}[(-1)^{s+1} p(n-\frac{3s^2+s}{2}) ])$

只需要初始化$p(0) = 1$

至此就已经可以编码了

有的地方会采取 $t=-s$的带入第一个求和，把上面两个求和合并成一个,也就是 s取非零的正负整数

$p(n) = \sum_{s\neq 0,s\in Z,\frac{3s^2+s}{2}\le n} [(-1)^{s+1}p(n-\frac{3s^2+s}{2})] $

然后就是五边形数相关的东西了,在这道题看来不是必须的

代码

fn main() {
    let mut ans = Vec::new();
    let modnum:i32 = 1_000_000; // 注意usize范围很小 大坑
    ans.push(1);
    let mut i = 1;
    loop{
        let mut ansi:i32 = 0;

        let mut s:i32 = 1;
        let mut sres = ( s*(3*s-1)/2 ) as usize;
        while sres <= i{
            ansi += ((s%2)*2-1) * ans[i - sres];
            ansi %= modnum;
            s+=1;
            sres = ( s*(3*s-1)/2 ) as usize;
        }

        s = 1;
        sres = (s*(3*s+1)/2) as usize;
        while sres <= i{
            ansi += ((s%2)*2-1) * ans[i - sres];
            ansi %= modnum;
            s+=1;
            sres = (s*(3*s+1)/2) as usize;
        }
        ansi+=modnum;
        ansi%=modnum;

        ans.push(ansi);
        println!("{}=>{}",i,ansi);
        if ansi == 0 {
            println!("ans :{}",i);
            return ;
        }
        i+=1;
    }
}

秒出

不过看结果和我第一遍卡爆内存的值一比，我想是不是多点内存也是暴力可行的XD。

本题没有用到的 euler’s theorem

互异分拆方案数=拆分成全是奇数个的方案数 ,$(1+x+x^2+…)(1+x^3+x^6+…)(1+x^5+x^{10}+…)…$

$\prod_{n=1}^{\infty} \sum_{k=0}^{\infty} x^{k*(2n-1)}$

$= \prod_{n=1}^{\infty} \frac{1}{1-x^{2n-1}}$

$= \prod_{n=1}^{\infty} \frac{1}{1-x^{2n-1}} \frac{\prod_{n=1}^{\infty} 1-x^{2n}}{\prod_{n=1}^{\infty} 1-x^{2n}}$

$= \frac{1}{\prod_{n=1}^{\infty} 1-x^{2n-1}}\frac{\prod_{n=1}^{\infty} (1-x^n)(1+x^n)}{\prod_{n=1}^{\infty} 1-x^{2n}}$

$ = \prod_{n=1}^{\infty} (1+x^n)$ (这。。。无限项这样拆分相消科学吗

就有了wikipedia上面那个公式 $\sum _{n=0}^{\infty }q(n)x^{n}=\prod _{k=1}^{\infty }(1+x^{k})=\prod _{k=1}^{\infty }{\frac {1}{1-x^{2k-1}}}.$

参考

https://en.wikipedia.org/wiki/Partition_(number_theory)

https://en.wikipedia.org/wiki/Pentagonal_number_theorem

题目

大意

解法

代码

参考

题目

大意