彩珠手串的配色计数

TSC999

又推出了S(m,7), S(m,8)的公式。
`\D S(m,7)=\frac{m(m^6+7m^3+6)}{14}`
`\D S(m,8)=\frac{m(m^7+4m^4+5m^3+2m+4)}{16}`

hujunhua

S(2,n)=A000029(n)
S(3,n)=A027671(n)
S(4,n)=A032275(n)
S(5,n)=A032276(n)
S(6,n)=A056341(n)
S(7,n)=A032276(n)

8色及以上的手链或手串数OEIS尚未收录，楼主若搞出结果来，可以申报上去。

S(m,9)=A060561(n)
S(m,10)=A060562(n)

kastin

用伯恩赛德引理（Burnside's lemma，波利亚计数定理就是用它推导得来的）就能得到如下结论：
对于 `m` 个颜色可重复地选取 `n` 个排成一圈，圆排列数为$$N(m,n)=\frac{1}{n}\sum_{i=1}^n m^{(n,i)}$$环排列数为
$$M(m,n)=N(m,n)/2+\begin{cases}\D\frac{m^{k+1}}{2},&n=2k+1\\
\D\frac{(m+1)m^k}{4},&n=2k
\end{cases}$$
经验证，符合1楼所给的经验公式。

hujunhua

对于`n`是一个奇素数` p` 的情况，我暗暗猜测的蒙对了。

`\D S(m,p)=\frac{m(m^{\frac{p-1}2}+p-1)(m^{\frac{p-1}2}+1)}{2p}`

hujunhua

对于其它一些特殊形式的珠数，也能得到显式解。
`\D S(m,2p)=\frac{m\left(m^{2p-1}+p\cdot m^p+(p+1)m^{p-1}+(p-1)(m+1)\right)}{4p}`
`\D S(m,2^t)=。。。。`

mathe

比如根据9#公式
$N(m,8)=1/8(m^8+m^4+2m^2+4m)$
$S(m,8)=1/16m(m+1)(m^6-m^5+m^4+3m^3+2m^2-2m+4)$

TSC999

哎哟，我原先还以为是个比较简单的问题，没想到踏进了 “波利亚” 的雷区，大大超出了本人的能力。

kastin

将二项式定理扩充至多项，引入`r`-重组合数 `\mathrm C_n^{n_1,n_2,\cdots,n_r}`，定义如下$$(x_1+x_2+\cdots+x_r)^n=\sum_{n_1+n_2+\cdots+n_r=n\\n_i\in \NN,\;i=1,2,\cdots,r}\mathrm C_n^{n_1,n_2,\cdots,n_r}x_1^{n_1}x_2^{n_2}\cdots x_r^{n_r}\tag{a}$$其中求和符号下面表示的是遍历所有可能满足条件的取值，`\NN` 为非负整数集。
利用二项式展开，有$$\mathrm C_n^{n_1,n_2,\cdots,n_r}=\mathrm C_n^{n_1}\mathrm C_{n-n_1}^{n_2}\mathrm C_{n-n_1-n_2}^{n_3}\cdots \mathrm C_{n_r}^{n_r}=\frac{n!}{n_1!n_2!\cdots n_r!}\tag{b}$$式 `(a)` 本身就是`r`-重集 `S(r)=\{n_1\*e_1,n_2\*e_2,\dots,n_r\*e_r\},\;(n_1+\cdots+n_r=n)` 的线排列的母函数，因此式 `(b)` 即为 `S(r)` 的线排列数 。
另外，令 `x_1=x_2=\cdots=x_r=1` 可得 $$r^n=\sum_{n_1+n_2+\cdots+n_r=n\\n_i\in \NN,\;i=1,2,\cdots,r}\mathrm C_n^{n_1,n_2,\cdots,n_r}=\sum_{n_1+n_2+\cdots+n_r=n\\n_i\in \NN,\;i=1,2,\cdots,r}\frac{n!}{n_1!n_2!\cdots n_r!}\tag{c1}$$这个结果很容易理解，左边表示 `n` 位数的每一位可以有 `r` 种选择，共 `r^n` 种选择，右端表示遍历所有可能重复的选择并求和。

通过3楼和9楼mathe的提醒，我们可以得到一个组合恒等式$$\sum_{d|n}\varphi(d)m^{\frac{n}{d}}=\sum_{n_1+n_2+\cdots+n_m=n\\n_i\in \NN,\;i=1,2,\cdots,m}\;\sum_{q|(n_1,n_2,\cdots,n_m)}\varphi(q)\frac{\frac{n}{q}!}{\frac{n_1}{q}!\frac{n_2}{q}!\cdots \frac{n_m}{q}!}\tag{c2}$$注意，上式右端中分子分母中的分式没有用圆括号括起来是因为避免与数论中的雅可比符号产生混淆（同样，3楼中对应位置本来也应该不加括号的）；还有右端之所以用 `q` 而不是用 `d` 表示约数，是为了避免相同符号在等式两边却有不同意义。

我感兴趣的是 `(\mathrm c1)` 和 `(\mathrm c2)` 之间的内在联系。易知，如果 `q` 遍历的值域就是 `n` 的所有因子（即和 `d` 的值域是一样的），那么这二式是等价的（注意到，二重求和符号可交换）。即问不定方程 `n_1+n_2+\cdots+n_m=n\;n_i\in \NN,\;(i=1,2,\cdots,m)` 的每组解的公约数构成的并集是否就是 `n` 的所有约数？

TSC999

根据 kastin 老师的公式，我用 mathematica 写了一段小程序，是当手串珠子数目 n 固定之后，配色方案数目 S(m, n) 的计算公式。
程序计算结果正确，但是不知为什么有警告提示。不知程序哪里不太好而导致，如何才能消除这个警告？
程序如下：

1. Clear[s]; Array[s, {30, 30}]; (*定义二维数组 s(16,16) *)
   
2. For[q = 1, q <= 30, q++,     (* 计算二维数组 s 的值 *)
   
3. For[p = 1, p <= 30, p++,
   
4.   If[EvenQ[q], s[p, q] = 1/(2 q) \!\(
   
5. \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(q\)]
   
6. \*SuperscriptBox[\(p\), \(GCD[q, i]\)]\) + ((p + 1) p^(q/2))/4,
   
7.     s[p, q] = 1/(2 q) \!\(
   
8. \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(q\)]
   
9. \*SuperscriptBox[\(p\), \(GCD[q, i]\)]\) + p^((q + 1)/2)/2];]]
   
10. n =.; j = 30;
    
11. Array[x, 30];             (* 解方程：用数组 x[i] 作为未知数 *)
    
12.   For[n = 1, n <= j, n++,(* 算 1 至 30 珠 *)
    
13. Clear[x]; a = Table[x[k], {k, 1, j}] /.
    
14.    Solve[Table[\!\(
    
15. \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(n\)]\(x[i]
    
16. \*SuperscriptBox[\(k\), \(i\)]\)\) == s[k, n], {k, 1, j}],
    
17.     Table[x[k], {k, 1, j}]];
    
18. b = a[[1]];
    
19. m =.; ss = \!\(
    
20. \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(n\)]\(b[\([\)\(i\)\(]\)]*
    
21. \*SuperscriptBox[\(m\), \(n + 1 - i\)]\)\);
    
22. Print["s(", n, ")=", Factor[Simplify[ss]]]]

复制代码

运行结果如下：

补充内容 (2021-7-11 18:30):
15# 这个程序的计算结果也有误。正确的结果见下面 16# 楼的程序给出的结果。

kastin

TSC999 发表于 2016-11-3 20:09
根据 kastin 老师的公式，我用 mathematica 写了一段小程序，是当手串珠子数目 n 固定之后，配色方案数目 S ...

Mathematica是函数式程序语言，学会函数式的编程风格，问题将会简单很多。

1. S[m_, n_?OddQ] :=
   
2.   1/2 Sum[EulerPhi[n/d] m^d, {d, Divisors[n]}]/n + m^((n + 1)/2)/2;
   
3. S[m_, n_?EvenQ] :=
   
4.   1/2 Sum[EulerPhi[n/d] m^d, {d, Divisors[n]}]/n + (m + 1) m^(n/2)/4;
   
5. Table[Defer[S][m, n] == S[m, n] // Factor, {n, 1, 6}] // TableForm

复制代码

以下是输出结果（必须在“编辑”>“偏好设置">”计算“>"新输出单元格式"，选择TranditionalForm)

输出结果