简介

hujunhua于2010年4月提出一个问题:
{1, 2, 3, …, n}的总体平均值a=(n+1)/2. 如果它的某个真子集的均值也是a, 就称为它的一个均衡样本。
一个均衡样本若包含更小的均衡样本，即为可约均衡样本，否则为不可约均衡样本。
问：{1, 2, 3, …, n}有多少个

1. k元均衡样本
2. k元不可约均衡样本
3. 不可约均衡样本

4. 极大划分（将{1, 2, 3, …, n}划分成若干个不可约均衡样本）

计算数据

hujunhua自己给出了第3）和第4）问的n=1~12的数据：

并且找出
3元序列为A000982
4元序列为A001973
6元序列为A001977

递推数列

wayne借助软件找出对于三元均衡样本, n应该为奇数, 样本个数为 $\lceil \frac{(n-1)^2}8 \rceil$.
并且随后给出更多的递推式结果:
(注释：K为奇数时，偶数n的样本个数为0，所以下面的3，5元反映的是去掉所有的偶数项的0之后的新序列)
三元: $a_n=2a_{n-1}-2a_{n-3}+a_{n-4}$
初始值{0, 1, 2, 5}
四元：$a_n=2a_{n-1}-a_{n-3}-a_{n-4}+2a_{n-6}-a_{n-7}$
初始值{0, 0, 0, 1, 1, 3, 5}
五元：$a_n=2a_{n-1}-a_{n-3}-2a_{n-5}+2a_{n-6}+a_{n-8}-2a_{n-10}+a_{n-11}$
初始值{0, 0, 1, 3, 12, 32, 73, 141, 252, 414, 649}
六元：$a_n=a_{n-1}+2a_{n-2}-a_{n-3}-a_{n-4}-a_{n-5}-a_{n-6}+2a_{n-8}+2a_{n-9}-a_{n-11}-a_{n-12}-a_{n-13}-a_{n-14}+2a_{n-15}+a_{n-16}-a_{n-17}$
初始值{0, 0, 0, 0, 0, 1, 1, 4, 8, 18, 32, 58, 94, 151, 227, 338, 480}

特征方程

hujunhua根据wayne的结果，给出了对应的特征多项式:
k=3时，$x^4 - 2 x^3 + 2 x - 1=(x-1)^3(x+1)=(x-1)^2(x^2-1)$
k=4时，$x^7 - 2 x^6 + x^4 + x^3 - 2 x + 1=(x-1)^4(x+1)(x^2+x+1)=(x-1)^2(x^2-1)(x^3-1)$
k=5时，$x^{11} - 2 x^{10} + x^8 + 2 x^6 - 2 x^5 - x^3 + 2 x - 1=(x-1)^5(x+1)^2(x^2+1)(x^2+x+1)=(x-1)^2(x^2-1)(x^3-1)(x^4-1)$
这些特征多项式的不可约因式都是分圆多项式，并且刚好可凑成1~(k-1)阶分圆方程之积，仅有因式(x-1)多了一重。
多么美妙的规律啊，他满怀期望地分解k=6的特征多项式，
$1 - x - 2 x^2 + x^3 + x^4 + x^5 + x^6 - 2 x^8 - 2 x^9 + x^{11}+ x^{12} + x^{13} + x^{14}-2x^{15} - x^{16} + x^{17}$
$=(x-1 )^6 (1 + x)^3 (1 + x^2) (1 + x + x^2) (1 + x + x^2 + x^3 + x^4)$
$=(x-1)(x^2-1)^2(x^3-1)(x^4-1)(x^5-1)$
结果这回多了一重的是$(x^2-1)$. 虽然仍然只有分圆多项式，仍然能够凑成1~(k-1)阶分圆方程之积，但是可惜重数已经没有规律了。

hujunhua进一步猜测：
记$Y_k(x)=\prod_{i=1}^{k-1}(x^i-1)$，k元均衡样本数序列的特征多项式记为$S_k(x)$。我们的合情猜想是：$Y_k(x)|S_k(x)$.
$Y_{7}(x)=1-x-x^2+x^{5}+2x^{7 }-x^{9}-x^{10}-x^{11}-x^{12}+2x^{14}+x^{16}-x^{19}-x^{20}+x^{21}$
构造$b_{n}=a_{n}-a_{1+n}-a_{2+n}+a_{5+n}+2a_{7+n}-a_{9+ n}- a_{10+n}-a_{11+n}-a_{12+n}+2a_{14+n}+ a_{16+n}-a_{19+n}-a_{20+n}+a_{21+n}$
b{n}满足的特征多项式就是剩余的因式。
得出$S{7}(x)=(x-1)Y{7}(x)$. 也许k=奇数时，都是$S{k}(x)=(x-1)Y_{k}(x)$.

mathe随后提供了一段Pari/gp计算代码

S(K,n)=
{
    local(V,L,s,R);
    L=K*n;
    V=matrix(K,L);
    R=vector(n);
    V[1,1]=1;
    R[1]=0;
    for(u=2,2*n-1,
       for(h=1,K-1,
          for(v=1,L,
            s=V[K+1-h,v];
            if(v>u,s=s+V[K-h,v-u]);
            V[K+1-h,v]=s
          )
       );
       V[1,u]=1;
       if(bitand(u,1)==1,
            s=(u+1)/2;
            R[s]=V[K,K*s]
       )
    );
    R
}

wayne利用mathe的代码验证了hujunhua的k=奇数时，都是$S_{k}(x)=(x-1)Y_{k}(x)$的猜想到15元都成立，但是偶数阶情况的规律不匹配。
他给出了部分偶数情况的特征方程:
6元： $(x^2-1)Y_6(x)$
8元：$\frac{-1+x}{1-x+x^2}Y_8(x)$
10元：$(x^2-1)Y_10(x)$
12元：$\frac{-1+x}{1+(-1+x) x (1+x^2)}Y_12(x)$
mathe猜测偶数阶特征方程都是$(x^2-1)Y_K(x)$的因子。

母函数

mathe发现,
根据A001973和A001977, 如果k是偶数，那么母函数应该是(包含偶数项时) $\prod_{s=1}^k\frac{x^s-x^n}{1-x^s}$.
并给出了计算母函数的pari/gp代码:

RR(K,n)=
{
   local(f,m,p,s,t);
   s=K+1;t=(K*(K+1))/2;
   f=1/(1-z);
   for(h=1,K,
     f=f/(1-x^h*z)
   );
   m=K*n;
   p=subst(subst(f,x,x+x*O(x^m)),z,z+z*O(z^m));
   for(h=1,n, if(2*h<K+1,print1(0","),print1(polcoeff(polcoeff(p,K*h-t),2*h-s)",")))
}

然后mathe假设母函数为$\frac{T_K(x)}{(x^2-1)Y_K(x)}$,
利用代码

S(K,n)=
{
    local(V,L,s,R);
    L=K*n;
    V=matrix(K,L);
    R=vector(n);
    V[1,1]=1;
    R[1]=0;
    for(u=2,2*n-1,
       for(h=1,K-1,
          for(v=1,L,
            s=V[K+1-h,v];
            if(v>u,s=s+V[K-h,v-u]);
            V[K+1-h,v]=s
          )
       );
       V[1,u]=1;
       if(bitand(u,1)==1,
            s=(u+1)/2;
            R[s]=V[K,K*s]
       )
    );
    R
}

EstPol(V,K)=
{
    local(len,p,r);
    len = length(V);
    p=V[1];
    for(u=2,len,
       p=p+V[u]*x^(u-1)
    );
    p=p*(1-x^2);
    for(u=1,K-1,
       p=p*(1-x^u)
    );
    r=polcoeff(p,0);
    for(u=1,len-1,
      r=r+polcoeff(p,u)*x^u
    );
    r
}
MS(k,n)=EstPol(S(k,n),k)

来计算$T_k(x)$.

wayne总结:
序列{1,2,3,...,n}的K元均衡样本的个数a(n)是$\prod_{s=1}^k\frac{x^s-x^n}{1-x^s}$的展开式的x 的 $\lfloor \frac{(n+1)k}2\rfloor$次幂的系数。
大一统的 Mathematica函数是：

f[kk_, nn_] := Module[{k = kk, tmp, n = nn}, tmp = PadLeft[Table[SeriesCoefficient[Series[Product[(x^s - x^(m + 1))/(1 - x^s), {s, 1, k}], {x, 0, m k}], Floor[(m + 1) k/2]], {m, k, n}], n]; If[EvenQ[k], tmp, tmp[[1 ;; -1 ;; 2]]]]

新的数列

多年以后，hujunhua重新来分析第三问:
他编写了一个M10小程序计算过第3）问，只计算到了前28项：
{1, 1, 2, 2, 3, 3, 6, 6, 11, 13, 20, 26, 41, 55, 74, 116, 141, 241, 258, 472, 473, 873, 828, 1576, 1447, 2821, 2472, 5072}
在OEIS上未搜到上述有限序列，应该是个新数列。
将奇数项和偶数项分开，得到两个单调递增数列：
{1, 2, 3, 6, 11, 20, 41, 74, 141, 258, 473, 828, 1447, 2472, ...}
{1, 2, 3, 6, 13, 26, 55, 116, 241, 472, 873,1576, 2821, 5072, ...}

为了程序的判断条件简明，代码中将自然数的前段 Range[n] 变换为 N(n)=Range[-n+1, n-1, 2], 例如
N(6)={-5, -3, -1, 1, 3, 5}
N(7)={-6, -4, -2, 0, 2, 4, 6}【程序中将0去掉了】
为的是使和与平均数皆为 0. 于是均衡样本即零和子集，不可约均衡样本即不含零和真子集的样本。
这个变换显示了长度为奇数的序列与长度为偶数的显然区别，故将两者分开处理是有道理的。
N(Odd)可以将其中的数除以2，下表为3-11的奇数按此处理的数据展示。可以与1#和4#的数据进行一下对照。

n(n+2)的不可约均衡样本可分为4种：1、N(n)的, 2、{-n-1,n+1}, 3、包含 -n-1而不包含 n+1 的, 4、包含 n+1 而不包含 -n-1 的.
其中后两种是负对称的，只需要计算二者之一。
前面说了，N(n+2)的不可约均衡样本可分为4种：1、N(n)的, 2、{-n-1,n+1}, 3、包含 -n-1而不包含 n+1 的, 4、包含 n+1 而不包含 -n-1 的.
记第4类的计数为d(n+2), 4类之和计数为a(n+2)，
则a(n+2)=a(n)+1+2d(n+2)=a(n-2)+2+2d(n)+2d(n+2)=...=a(5)+(n-3)/2+2d(7)+2d(9)+...+2d(n+2)=(n+3)/2+2d(7)+2d(9)+...+2d(n+2).
我们的程序就是搜索计数d(k), 然后累加得到各a(n).

1. 生成一个长为(n-1)/2的对称三进数位表, 如{0,1,-1,-1,0,0,...,1}. 函数为Tuples[{-1, 0, 1}, (n-1)/2]
2. 将上述位表与序列{2,4,6,...,n-1}对位相乘，滤掉零元素。得到对序列的选取和赋号，然后将n+1加入得表最后。
3. 滤取得表集中的零和列表（第一步中的函数得到的其实是所有的位表集）。

4. 构造一个检验函数， 检验所得表列是否包含零和真子集。用此函数滤掉包含零和真子集者。

   既约均衡样本数

   n	既约均衡样本（省略0）
   13	41	{±1},{±2},{±3},{±4},{±5},{±6}
   		±{1,2,-3},±{1,3,-4},±{1,4,-5},±{1,5,-6},
   		±{2,3,-5},±{2,4,-6},
   		±{1,-2,-3,4},±{1,2,3,-6},±{1,-2,-4,5},±{1,-2,-5,6},
   		±{1,-3,-4,6}±{2,-3,-4,5},±{2,-3,-5,6},±{3,-4,-5,6}
   		±{1,-2,3,4,-6},±{1,2,-4,-5,6},±{2,3,-4,5,-6}