特殊数列

最后更新于：2023年2月24日 下午

在知乎上看到一个有趣的问题: 如何证明这个数列无界？ 在此记录一下简单的做法，然后把这篇博客记录以后遇到的一些特殊数列。

证明：当 \(x\) 是无理数时，\(f_n(x) = \sum_{i=1}^n (-1)^{\lfloor ix \rfloor}\) 无界

\(\lfloor x \rfloor\) 表示 \(x\) 的整数部分， \(0 \leq \lbrace x \rbrace \doteq x- \lfloor x \rfloor<1\) 表示 \(x\) 的小数部分

为了更顺畅的介绍这个数列，我们先来点前戏，不然插入的时候就不那么顺滑了 0.0 不要笑，我没有开车 0.0

连分数

它和 Farey 序列，无理数的有理逼近，密切相关，还是由于求 Pell 方程的快速方法

我们将一个整数序列 \(a_0,a_1, \cdots, a_n\) 构成的分数叫做连分数，记作 \[ [a_0,a_1,\cdots, a_n] = a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \cdots}} \] 其中 \(a_0\) 是整数， \(a_i (i>0)\) 是正整数。

我们定义:

\(p_0 = a_0,p_1 = a_1a_0 +1,p_n = a_np_{n-1} + p_{n-2} (n>1)\)

\(q_0 = 1,q_1 = a_1,q_n=a_nq_{n-1}+q_{n-2} (n>1)\)

那么数学归纳法易证： \[ [a_0,a_1,\cdots, a_n] = \frac{p_n}{q_n} = x_0 \] 并且还有

  • \([a_0,a_1,\cdots, a_n] = [a_0,a_1,\cdots, a_{n-1} + \frac{1}{a_n}]\)
  • \(p_{n+1}q_n - q_{n+1}p_n = (-1)^n\)
  • 从而 \((p_n,p_{n+1}) = (q_n,q_{n+1}) = (p_n,q_n) = 1\)
  • \(x_{n+2} - x_n = \frac{(-1)^n a_{n+2}}{q_n q_{n+2}}\)
  • 从而 \(x_0 < x_2< x_4 < \cdots < x_5 < x_3 < x_1\)

对任意一个数 \(x\) 构造连分数

  1. 初始 \(i = 0\)
  2. \(a_i = \lfloor x \rfloor\) 取 \(x\) 的整数部分。
  3. \(x - a_i == 0\) 结束
  4. 取 \(x = \frac{1}{x-a_i}\)
  5. ++\(i\) 回到步骤 1

从构造中，我们看出：连分数长度有限当且仅当 \(x\) 是有理数，此时，最后一项得到的连分数就是 \(x\)。

当 \(x\) 是无理数时， \(x_{2n} < x < x_{2n+1}\) ，又因为 \[ x_{2n+1} - x_{2n} = \frac{1}{q_{2n}q_{2n+1}} \leq \frac{1}{2n(2n+1)} \] 且 \(x_{2n}\) 单调递增，\(x_{2n+1}\) 单调递减。所以 \(\lim _{n \to \infty} x_n = x\)。并且 \[ |x - \frac{p_n}{q_n} | < \frac{1}{q_n ^2} \]

补充定义

\[ a_n' \doteq a_n'(x) \doteq [a_n, a_{n+1}, a_{n+2},\cdots] \]

\[ p_{n+1}' \doteq a_{n+1}' p_n + p_{n-1} = \frac{p_{n+2}'}{a_{n+2}'} \\ q_{n+1}' \doteq a_{n+1}' q_n + q_{n-1} = \frac{q_{n+2}'}{a_{n+2}'} \]

显然，\(x = \frac{p_n'}{q_n'} = \frac{p_n}{q_n} + \frac{(-1)^n}{q_n q_{n+1}'}\)

唯一分解

对于任意给定 \(x = [a_0,a_1,\cdots,a_n,\cdots]\)，对任意 \(N\)，有如下唯一分解： \[ N = \sum c_i q_i = \sum_{i=0} ^{m(N,x)} c_i(N,x)q_i(x) \] 其中 \(0 \leq c_i \leq a_{i+1}, c_m(N,x) >0\) 且对任意 \(0 \leq j \leq m(N,x)\)，\(\sum_{i=0} ^ j c_i q_i < q_{j+1}\) 。

$ f_n(x) = _{i=1}^n (-1)^{ix } $ 何时有界

由于 $f_n(x+2) = f_n(x) $， 所以我们仅考虑 $x $ 时，\(i\frac{p}{2} \mod q\) 跑遍 \([0,q-1]\)，即 \[ f_q = \sum_{i=1}^q (-1)^{ \lfloor \frac{2( i\frac{p}{2} \mod q)}{q} \rfloor} =\sum_{i=0} ^{q-1} (-1)^{\lfloor \frac{2i}{q} \rfloor} = \frac{q+1}{2} - \frac{q-1}{2} = 1 \] 从而 \(f_{2q} = 2f_q = 2\) 推出 \(f_n\) 无界。

对任意无理数 \(x\)，我们证明 \(g_n(x) = f_n(2x)\) 无界，从而 \(f_n(x)\) 无界

目前没理解论文中的做法，太繁琐了

显然，我们只需考虑 \(0< x <1\)

我们考虑 \(x\) 的连分数 \([0,a_1,a_2,\cdots,a_n,\cdots]\)，由于 \((p_n,q_n)=(q_n,q_{n+1})=1\) (目前不知道有什么简单的操作)

这个问题可以推广到更一般的情形：单栏阅读 双栏阅读

里面定理一挺有意思的，虽然长但是很清晰。

另一个相关问题：https://www.zhihu.com/question/392014769

是否存在 \(N\) 使得， \(|\sin n|>\frac{1}{n}\) 对所有 \(n>N\) 成立。(cuzperf 和 寨森 Lambda-CDM 建议修改)

存在无穷多个正整数 \(x_n\) 使得 \(|\sin x_n| < \frac{2 \pi}{x_n}\)

引理： 对任意无理数 \(a\) , 存在无穷多个 \(x_n\)，使得 \(\min(\{a x_n\},\{-a x_n\}) < \frac{2}{x_n}\)

我们考虑 \(a\) 的连分数 \([a_0,a_1,\cdots,a_n,\cdots]\), 我们知道 \(|a - \frac{p_n}{q_n} | < \frac{1}{q_n ^2}\)，由于 \((p_n,q_n) = 1\), 所以存在 \(u_n,v_n\) 使得 \(p_n u_n + q_n v_n = 1\), 此时 \(p_n(u_n+tq_n) + q_n(v_n - tp_n) = 1\)，所以我们不妨假设 \(q_n < u_n < 2q_n\)。所以 \[ \{ a u_n \} = \{ (a-\frac{p_n}{q_n})u_n + \frac{p_n}{q_n}u_n \} = \{ (a-\frac{p_n}{q_n})u_n + \frac{1}{q_n} \} \] 我们仅考虑 \(n\) 为奇数数的情况，此时 \(-\frac{1}{q_n^2} < a - \frac{p_n}{q_n} < 0\)，所以 \(-\frac{1}{q_n} <(a-\frac{p_n}{q_n})u_n + \frac{1}{q_n} <\frac{1}{q_n}\)。而 \(x_n = u_{2n+1} > q_{2n+1}\) 即为所求。

实际上 \(n\) 为偶数时，也可以取做，只是此时 \(u_n\) 要换成 \(3q_n - u_n\)

取引理中 \(a = \frac{1}{\pi}\)，由于 \(|\sin(x)|\) 是周期为 \(\pi\) 的偶函数。所以 \(|\sin(x_n)| = |\sin(\{ax_n\}\pi)| = |\sin(\{-ax_n\}\pi)| < \frac{2\pi}{x_n}\)

受 寨森 Lambda-CDM 启发，可以不用引理直接证明：

存在无穷多个正整数 \(m\)，使得 \(|\sin m| < \frac{\pi}{m}\)

Proof：考虑 \(\frac{1}{\pi}\) 的连分数 \([a_0,a_1,\cdots,a_n,\cdots]\)，我们有 \(|\frac{1}{\pi} - \frac{p_n}{q_n} | < \frac{1}{q_n ^2}\)，所以 \[ |\sin q_n| = |\sin[(\frac{q_n}{\pi} -p_n) \pi + p_n\pi ]| = |\sin[(|\frac{q_n}{\pi} -p_n|) \pi]| < \frac{\pi}{q_n} \]

这个问题还关联这 Irrationality Measure，菲尔兹奖级别的工作！

Irrationality Measure

对于给定是实数 \(x\), 定义 \[ \begin{aligned} \mu(x) \doteq \inf_{u \in R} u \\ R = \{ u\mid \exists \text{ infty } (p,q)\; s.t.\quad 0< |x -\frac{p}{q}| < \frac{1}{q^u} \} \end{aligned} \] 显然考虑连分数，我们知道 当 \(x\) 是有理数时，\(u(x) = 1\)，无理数时，\(\mu(x) \geq 2\)。

Roth 证明，当 \(x\) 是代数数时 \(\mu(x) = 2\)，因此获得 Field 奖。

\(\mu(L) = \infty\)，其中 \(L = \sum_{n=1} ^{\infty} 10^{-n!}\) 为刘维尔（Joseph Liouville）数

我们下面证明：当 \(u>2\) 时，\(A = \{ x \mid 0< |x -\frac{p}{q}| < \frac{1}{q^u}\}\)，则 \(A\) 的外（lebesgue）测度 \(m^{\star}(A) = 0\)，从而（lebesgue）测度 \(m(A)=0\)

由于 当 \(1 \leq p \leq q\) 时， \((\frac{p}{q} - \frac{1}{q^u}, \frac{p}{q} + \frac{1}{q^u})\) 只有可数个，记作 \(I_1, I_2,\cdots I_n,\cdots\) ，按照定义，\(A \cap (0,1) \subseteq \overline\lim I_n\), 另一方面 \(\sum_{i=1}^{\infty} I_n = \sum_{q=1}^{\infty} \frac{q+1}{q^u} < +\infty\)，从而 \(m^{\star}(A \cup (0,1)) = 0\)，同理可证了\(m^{\star}(A \cup (n,n+1))=0, n \in \mathbb{Z}\)，从而 \[ 0 \leq m^{\star}(A) \leq \sum_{n \in \mathbb{Z}} m^{\star}(A \cup (n,n+1)) =0 \] 从而\(m^{\star}(A) = m(A) = 0\)。其中 \(m^{\star}(I) = \inf \{u|u=\sum_{k=1}^{\infty} |I_k|,\; \cup_{k=1}^ {\infty} I_k \supset m ,I_k \text{ 是开矩形 } \}\)

sagemath 数值测试(Jupyter Notebook 上运行)

1
2
3
4
5
6
7
8
9
x = continued_fraction(1/pi).convergents()
for i in range(1,32):
    n = x[i].denominator()
    print(N(abs(sin(n))-pi/n),"\t= ",abs(sin(n))-pi/n ) # 后来精度就崩了

x = continued_fraction(1/pi).convergents()
for i in range(1,32):
    n = x[i].denominator()
    print(i,N(abs(sin(n))*n),"\t= ",abs(sin(n))*n) # 后来精度就崩了

\(1 = \sum_{s \in S} \frac{1}{s^2-1}\), 其中 \(S = \{a^m | a > 1, m>1, a,m \in \mathbb{N} \}\)

方法：On a Series of Goldbach and Euler

\[ \begin{aligned} \lim S(n) &= \sum_{a \in A} \sum_{m\geq 2} \frac{1}{a^m -1} \\ &= \sum_{a \in A} \sum_{m\geq 2} \frac{1}{a^m} + \sum_{a \in A} \sum_{m\geq 2} \frac{1}{a^m(a^m -1)} \\ & = \sum_{a \in A} \frac{1}{a(a-1)} + \sum_{b \in B} \frac{1}{b(b -1)} \\ & = \sum_{i=2}^n \frac{1}{n(n-1)} = 1 \end{aligned} \]

其中 \(B = \{n \in \mathbb{N} \mid n = k^m , k>1 \}\)，\(A= \{n \in \mathbb{N} \mid n >1 \} - B\)

\(f(n,m) = \sum_{i=0}^{i \leq n-im} C_{n-im}^i\)

洛谷题目，具体题号就不便说了。

数据范围 \(1 \leq n \leq 10^{18} , m \leq 100\) 。我们注意考察 \(f(n,m)\) 的意义！假设我们要上楼梯，每次只能上 1 步或者 \(m+1\) 步，那么 \(f(n,m)\) 就是方案数，因此，我们显然有

\[ f(n,m) = f(n-1,m)+f(n-m-1,m) \]

当然了我们不考虑意义直接用 \(f(n,m)-f(n-1,m)\) 也可以得到这个公式。

无法过题的 \(O(m \log m \log n)\) 代码，原因是没法找一个合适的基底

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#include<bits/stdc++.h>
using namespace std;
using LL=long long;
using BI=__int128;
const LL mod=1e9+7;
const int M =10004;
BI a[4*M],r[4*M],ans[2*M];
template<typename T,typename U>
T powmod(T x,U n,T p){
  T r(1);
  while(n){
    if(n&1) r=r*x%p;
    n>>=1;  x=x*x%p;
  }
  return r;
}
void bitreverse(BI *x,int len){
  for(int i=0,j=0;i!=len;++i){
    if(i>j) swap(x[i],x[j]);
    for(int l=len>>1;(j^=l)<l;l>>=1);
  }
}
const BI FM = BI(29)<<57|1,gg=3;
// the mod must NTT-friendly or (len+1)*mod^2 < FM
void ntt(BI *x,int len,bool isInverse=false){
  g = powmod(gg,(FM-1)/len,FM);
  if(isInverse){
    g=powmod(g,len-1,FM);
    BI invlen = powmod(BI(len),FM-2,FM);
    for(int i=0;i!=len;++i){
      x[i]=x[i]*invlen%FM;
    }
  }
  bitreverse(x,len);
  for(int half=1,step=2;half!=len;half<<=1,step<<=1){
    BI wn = powmod(g,len/step,FM);
    for(int i=0;i<len;i+=step){
      BI w(1);
      for(int j = i;j<i+half;++j){
        BI t=(w*x[j+half])%FM;
        x[j+half]=(FM-t+x[j])%FM;
        x[j]=(x[j]+t)%FM;
        w = w*wn%FM;
      }
    }
  }
}
void stand(BI *a, int n){
  for(int i=2*n;i>n;--i){
    a[i-1]=(a[i-1]+a[i])%mod;
    a[i-n-1]=(a[i-n-1]+a[i])%mod;
    a[i] = 0;
  }
}
void square(BI *a, int n){
  int len = 1<<(32-__builtin_clz(2*n+1));
  ntt(a,len);
  for(int i=0;i<len;++i){
    a[i]=a[i]*a[i]%FM;
  }
  ntt(a,len,1);
  stand(a,n);
}
void mul(BI *a,BI *b,int na,int nb){
  int len = 1<<(32-__builtin_clz((na+nb)+1));
  ntt(a,len);ntt(b,len);
  for(int i=0;i<len;++i){
    a[i] = a[i]*b[i]%FM;
  }
  ntt(b,len,1);ntt(a,len,1);
  stand(a,na);
}
void initC(int m){
  for(int i=0;i<=m;++i){
    ans[i] = 1;
  }
  for(int i=m+1;i<=2*m;++i){
    ans[i] = (ans[i-1]+ans[i-m-1])%mod;
  }
  memset(a,0,sizeof(a));
  memset(r,0,sizeof(r));
  r[0]=a[1]=1;
}
LL getans(LL n,int m){
  initC(m);
  while(n){
    if(n&1) mul(r,a,m,m);
    n>>=1; square(a,m);
  }
  LL ret = 0;
  for(int i=0;i<=m;++i){
    ret+=(r[i]*ans[i+m])%mod;
  }
  return ret%mod;
}

int main() {
  //freopen("in","r",stdin);
  //freopen("m","w",stdout);
  std::ios::sync_with_stdio(false);std::cin.tie(nullptr);
  LL n;
  int m;
  while(cin>>n>>m){
    if(m==0){
      cout<<powmod(2LL,n,mod)<<endl;
    }else{
      if(n<=m){
        cout<<1<<endl;
      }else{
        cout<<getans(n-m,m)<<endl;
      }
    }
  }
  return 0;
}

求 \(m\) 阶线性递推关系式第 \(n\) 项

参考：2013/03/02 线性递推关系和矩阵乘法

上面文献的意义：矩阵加速已经成为过去式了，多项式加速时代的到来，NTT 加速多项式的时代到来。

设 \(a_n = c_{m-1} a_{n-1} + c_{m-2} a_{n-2} + c_0 a_{n-m}\)，给定了初值 \(a_1,a_2,\cdots,a_m\) 的情况下，求 \(a_n\)

我们记 \[ A = \begin{pmatrix} c_1 & c_2 &\cdots& c_m \\ 1 & 0& & 0 \\ & &\vdots& \\ 0 & \cdots &1 & 0 \end{pmatrix} \] 并且显然 \[ \begin{pmatrix} a_n \\ a_{n-1} \\ \vdots \\ a_{n-m+1} \end{pmatrix} = A \begin{pmatrix} a_{n-1} \\ a_{n-2} \\ \vdots \\ a_{n-m} \end{pmatrix} = A^{n-m} \begin{pmatrix} a_{m} \\ a_{m-1} \\ \vdots \\ a_{1} \end{pmatrix} \] 我们就可以用矩阵乘法在 \(O(m^3 \log n)\) 解决这个问题了。

但线性递推关系式可以不用矩阵优化到 \(O(m^2 \log n)\) ：

显然 \(A\) 的特征多项式 \(f(x) = x^m - c_1 x^{m-1}- \cdots -c_m\)，所以 \(A^n = c_{m-1} A^{n-1} + \cdots c_0 I\)，也就是说 \(A^n\) 可以由 \(A^{m-1},A^{m-2},\cdots,I\) 线性表出，而 \(A^{n_1+n_2} = A^{n1} A^{n_2}\) 然后运算就是卷积运算，相当于多项式乘法！也就是说我们可以在 \(O(m^2 \log n)\) 时间复杂度，求出 \(A^n = b_0 I + b_1 A + \cdots b_{m-1} A^{m-1}\), 注意到 \[ A^i \begin{pmatrix} a_{m} \\ a_{m-1} \\ \vdots \\ a_{1} \end{pmatrix} = \begin{pmatrix} a_{i+m} \\ a_{i+m-1} \\ \vdots \\ a_{i+1} \end{pmatrix} \] 所以我们只需预处理出 \(a_1,\cdots,a_{2m-1}\) 这 \(2m-1\) 个数即可。 代码见我博客的模板

这也提供了一般 \(m\) 阶矩阵的 \(n\) 次方的一个 \(O(m^3+m^2\log n)\) 的算法！！！卧槽！我也太帅了吧。

最后如果递推关系中仅有常数个 \(c_i\) 不为 0，此时还能用 NTT（数论快速变换） 利用 多项式带模除法: Division with remainder 的 \(O(m \log m)\) 算法，上述算法可以优化到 \(O(m\log m \log n + m^2)\)，（暂时不知道如何去掉 \(m^2\)）但是写法太复杂。搞定！参考 这里，不要涉及矩阵

注意到一次乘法之后，会变成 \(I,A,\cdots A^{2m-1}\) 的线性组合，\(A^{m},\cdots A^{2m-1}\) 这 \(m\) 个要再用前 \(m\) 线性表出，由于 \(c_i\) 仅有常数个不为 0，可以 \(O(m)\) 复杂度把他们写成前 \(m\) 个的线性组合。

当 \(m>10^4\) 时，我们就有必要写带 NTT 的版本了（NTT 模板可以在我博客中找到）有需求的时候再写吧。

\(\sum_{i \equiv r \mod m} \binom{n}{i} \mod M\)

数据范围：\(1 \leq n \leq 10^{18},\; 2 \leq m \leq 2000,\; 0 \leq r < m,\; 10^8 < M < 10^9\)

记 \(w\) 满足 \(w^m = 1,w^n \neq 1, 0 < n < m\)，则答案为 \(\frac{1}{m}\sum_{i=0}^{m-1}F(w^i)\)，其中 \(F(x) = x^{m-r}(1+x)^n \mod x^{m}-1\)

假设 \(F(x) = \sum_{i=0}^{m-1} a_i x^i\)，则答案就是 \(a_0\)，即答案是 \(F(0)\)

这个跟上面一样本质是一样的，做带模的多项式运算。

参考：zscoder 的博客

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include<bits/stdc++.h>
using namespace std;
using LL = long long;

class Poly{
 public:
  const static int N = 2003;
  int d;
  LL a[N];
  Poly(int _d):d(_d){
    memset(a,0,sizeof(a));
  }
  auto& operator[](int i){
    return a[i];
  }
  Poly operator*(Poly &A){
    Poly R(d+A.d);
    for(int i=0;i<=d;++i){
      for(int j=0;j<=A.d;++j){
        R[i+j] += a[i]*A[j];
      }
    }
    return R;
  }
};
Poly mulmod(Poly &A,Poly &B,int m,LL mod){
  Poly R = A*B;
  for(int i=R.d;i>=m;--i){
    R[i-m] += R[i];
    R[i]=0;
  }
  R.d = m;
  for(int i=0;i<=m;++i){
    R[i]%=mod;
  }
  return R;
}
template<typename T,typename U>
T powmod(T x,U n,int m,LL mod){
  T r(0);
  r[0]=1;
  while(n){
    if(n&1) r=mulmod(r,x,m,mod);
    n>>=1;  x=mulmod(x,x,m,mod);
  }
  return r;
}

int main(){
  //freopen("in","r",stdin);
  std::ios::sync_with_stdio(false);std::cin.tie(nullptr);
  LL n,r,mod;
  int m;
  while(cin>>n>>m>>r>>mod){
    Poly A(m-1);
    A[0]=A[1]=1;
    Poly R = powmod(A,n,m,mod);
    cout<<(R[r%m]*m%mod)<<endl;
  }
  return 0;
}

组合数倒数交错和 \(\sum_{i=0}^n \frac{(-1)^i}{n \choose i}\)

考虑积分 \(F(n,m) = \int_0^1 x^m (1-x)^{n-m} dx\)，显然 \(F(n,m)=F(n,n-m),F(n,0) = \frac{1}{n+1}\) 且分部积分即可知道 \(F(n,m) = \frac{m}{n-m+1} F(n,m-1) = \frac{1}{(n+1){n \choose m}}\)，再等比数列求和就有

\[ \sum_{i=0}^n \frac{(-1)^i}{n \choose i} = (n+1) \int _0 ^1 \sum _{i=0} ^n (-x)^ i (1-x)^ {n-i} dx = (n+1) \int_0 ^1 (1-x)^{n+1} - (-x)^{n+1} dx = \frac{n+1}{n+2} (1+(-1)^n) \]

To be continued