系数为0,1的n次多项式的实根数

hbghlyj

给定正整数$n$,对于任意$S\subset\{0,\dots,n\}$,方程$$\sum_{i\in S} x^{i}=0$$的不同实根数的最大值$a_n$

$n$	1	2	3	4	5	6	7	8	9	10	11	12	13	14
$a_n$	1	2	2	2	2	2	3	3	3	4	4	4	4	4

$n=7$时

1. Max[Length@DeleteDuplicates@NSolve[Total[x^#]==0,x,Reals]&/@Subsets[Range[0,7]]]

Copy the Code

输出3
例如$x^7+x^4+x^2+x=0$

hbghlyj

去除$f(x)$的重根$${\rm rad}(f(x))\, =\, \frac{f(x)}{\gcd(f(x),f'(x))}$$
例如

1. syms x;
2. f=x^2-2*x+1;
3. [~,q]=polynomialReduce(f,diff(f));
4. q

Copy the Code

输出 x/2 - 1/2
用上面的写成Matlab代码:

1. syms x;
2. for i=2:11
3.     display([i-1 maxroots(i)]);
4. end
5. function max_len=maxroots(n)
6.     max_len = 0;
7.     combinations = dec2bin(0:2^n-1) - '0';
8.     for i = 2:2^n
9.         poly=poly2sym(combinations(i,:));
10.         [~,q]=polynomialReduce(poly,gcd(poly,diff(poly)));
11.         r=root(q);
12.         len=length(r(vpa(imag(r))==0));
13.         if len > max_len
14.             max_len = len;
15.         end
16.     end
17. end

Copy the Code

结果:

$n$	$a_n$
1	1
2	2
3	2
4	2
5	2
6	2
7	3
8	3
9	3
10	4

和1#算的结果一样. 对$n=11$计算时间太长就停止了. 有没有更高效的算法?

hbghlyj

参考MSE
使用Sturm algorithm改写一下:

1. for i=2:25
2.     display([i-1 maxroots(i)]);
3. end
5. function max_len=maxroots(n)
6.     max_len = 0;
7.     combinations = dec2bin(1:2^n-1) - '0';
8.     for i = 1:2^n-1
9.         c = combinations(i,:);
10.         if sum(c, 2) == 1
11.             continue;
12.         end
13.         len = distinct_real_roots(c);
14.         if len > max_len
15.             max_len = len;
16.         end
17.     end
18. end
19. function sign_var=distinct_real_roots(a)
20. S=Sturm();
21. S.build(Poly(a));
22. sign_var=0;
23. last_sign=0;
24. last_sign_parity=0;
25. parity=1;
26. for i=1:length(S.m_sturm)
27.     v=S.m_sturm{i}.m_coeffs(end);
28.     if v>0
29.         if last_sign<0
30.             sign_var = sign_var - 1;
31.         end
32.         if last_sign_parity == -parity
33.             sign_var = sign_var + 1;
34.         end
35.         last_sign = 1;
36.         last_sign_parity = parity;
37.     elseif v<0
38.         if last_sign>0
39.             sign_var = sign_var - 1;
40.         end
41.         if last_sign_parity == parity
42.             sign_var = sign_var + 1;
43.         end
44.         last_sign = -1;
45.         last_sign_parity = -parity;
46.     end
47.     parity = -parity;
48. end
49. end

Copy the Code

比2#的代码高效. 等到运行到15时间很长,我就把它停止了:

$n$	$a_n$
1	1
2	2
3	2
4	2
5	2
6	2
7	3
8	3
9	3
10	4
11	4
12	4
13	4
14	5
15	5

OEIS搜索结果有2个
oeis.org/search?q=1+2+2+2+2+2+3+3+3+4+4+4+4+5+5&go=Search

hbghlyj

hbghlyj 发表于 2023-5-5 01:42
    combinations = dec2bin(1:2^n-1) - '0';

$a_n$为$≤n$次0,1多项式的不同实根数最大值
这个代码计算$a_n$时包含了重复计算$a_1,\dots,a_{n-1}$
令$b_n$为$n$次0,1多项式的不同实根数最大值
那么$a_n=\max_{m≤n}b_m$
为计算$a_n$我们只需计算$b_n$
dec2bin(1:2:2^3-1)-'0'

ans =

     0     0     1
     0     1     1
     1     0     1
     1     1     1

dec2bin(1:2:2^4-1)-'0'

ans =

     0     0     0     1
     0     0     1     1
     0     1     0     1
     0     1     1     1
     1     0     0     1
     1     0     1     1
     1     1     0     1
     1     1     1     1
当$n>3$时($b_3=2$当$x^3+x^2$取最大值2),我们可丢弃被$x^2$整除(即第1,2个系数为0)的多项式$f(x)$,因为$f(x)$的不同实根数=多项式$f(x)\over x$的实根数.
dec2bin(2^(3-2)+1:2:2^3-1)-'0'

ans =

     0     1     1
     1     0     1
     1     1     1

dec2bin(2^(4-2)+1:2:2^4-1)-'0'

ans =

     0     1     0     1
     0     1     1     1
     1     0     0     1
     1     0     1     1
     1     1     0     1
     1     1     1     1
只改了这行,和i的范围,其它代码不变:

1. for i=4:20
2.     display([i-1 maxroots(i)]);
3. end
5. function max_len=maxroots(n)
6.     max_len = 0;
7.     combinations = dec2bin(2^(n-2)+1:2:2^n-1) - '0';
8.     for i = 1:3*2^(n-3)
9.         c = combinations(i,:);
10.         if sum(c, 2) == 1
11.             continue;
12.         end
13.         len = distinct_real_roots(c);
14.         if len > max_len
15.             max_len = len;
16.         end
17.     end
18. end
19. function sign_var=distinct_real_roots(a)
20. S=Sturm();
21. S.build(Poly(a));
22. sign_var=0;
23. last_sign=0;
24. last_sign_parity=0;
25. parity=1;
26. for i=1:length(S.m_sturm)
27.     v=S.m_sturm{i}.m_coeffs(end);
28.     if v>0
29.         if last_sign<0
30.             sign_var = sign_var - 1;
31.         end
32.         if last_sign_parity == -parity
33.             sign_var = sign_var + 1;
34.         end
35.         last_sign = 1;
36.         last_sign_parity = parity;
37.     elseif v<0
38.         if last_sign>0
39.             sign_var = sign_var - 1;
40.         end
41.         if last_sign_parity == parity
42.             sign_var = sign_var + 1;
43.         end
44.         last_sign = -1;
45.         last_sign_parity = -parity;
46.     end
47.     parity = -parity;
48. end
49. end

Copy the Code

输出

$n$	$b_n$
4	2
5	2
6	2
7	3
8	3
9	3
10	4
11	4
12	4
13	4
14	5
15	5
16	5

hbghlyj

OEIS draft

10:35	Robert Israel: Is there any reason to think this should be the same as A090735 or A090736, other than the coincidence of computed values? I don't see why they should be related.
11:03	Robert Israel: Can you give an example with 5 roots for n = 14, 15 or 16? My Maple program only gets 4.
12:45	Robert Israel: The first appearance of 5 roots I get is at n = 19, with x^19+x^18+x^16+x^11+x^9+x^6+x^5+x^4+x^2+x.
12:56	Robert Israel: So my program gets 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5. Can anyone confirm?