手串上穿有六颗珠子，求各珠子上的数字。

mathe

pari/gp计算33棵珠子的方法
? u=ffinit(2,5)
%82 = Mod(1, 2)*x^5 + Mod(1, 2)*x^4 + Mod(1, 2)*x^2 + Mod(1, 2)*x + Mod(1, 2)
? a=ffgen(u,a) //产生32阶有限域生成元，上面一行是对应生成多项式
%83 = a
? factor(x^3+a*x+a) //手工试验出一个3次本原多项式，这一步验证不可以因子分解
%84 =
[x^3 + a*x + a 1]
? factor(2^15-1) //因子分解$F_{r^3}$的阶
%85 =
[  7 1]
[ 31 1]
[151 1]
? Mod(x,x^3+a*x+a)^((2^15-1)/7) //依次计算Mod(x,x^3+a*x+a)的(2^15-1)除以7,31,151次不是1证明x^3+a*x+a是本原多项式
%86 = Mod((a^4 + a + 1)*x^2 + a^2*x + 1, x^3 + a*x + a)
? Mod(x,x^3+a*x+a)^((2^15-1)/31)
%87 = Mod(a, x^3 + a*x + a)
? Mod(x,x^3+a*x+a)^((2^15-1)/151)
%88 = Mod((a^2 + a + 1)*x^2 + (a^4 + a^2 + a)*x + a^3, x^3 + a*x + a)

?for(u=1,(2^15-1)/31, ;if(polcoeff(component(Mod(x,x^3+a*x+a)^u,2),2)==0,print(u))) //计算x^2系数为0的点的次数（也就是无穷远直线），得出
1
3
7
15
31
63
127
150
172
208
250
255
301
326
345
417
501
511
528
603
614
632
653
691
792
835
844
924
950
990
1003
1023
1057
这些数据相邻数的差就构成了一个33棵珠子的手链（其中最后一个数1057相当于0）

白新岭

“里面说明了如果n≡1(mod4)或者n≡2(mod4)而且n包含4k+3形式的素因子，那么不存在长度为n+1的手串”
以上是30楼结论，但它与实际情况不符，例如11无解，11-1=10，10≡2(mod4)但是它不含4k+3素数因子，也就是说它不符合无解的条件，可是无解；
例如13无解，13-1=12，12≡0(mod4)，12=4*3含有4k+3形式的素数因子，只符合条件之一(或为选择，算一种条件，而且n包含4k+3形式的素因子，算一种条件）。
从这两个例子看，那条件的关系都是“或”，对于自然数而言，如果按模4来说，只能把它分成4类，有两成已经排除，还有一个问题即其因子含4k+3，整除数中的因子只有3种情况，一个含素数因子4k+1，另一种含4k+3，第三种两种因子都含，这三种情况1，2的各占1/4,而后一种占2/4,所以又排除了1/4*3/4=3/16;   而含有4k+3因子，与n≡3(mod4)是一致的，所以只有1/16的自然数n加1后是有解的，这与实际情况不符。

白新岭

N=10，10-1=9=3*3,  9≡1(mod4)而且含有4k+3的素数因子3，但是它有解，所以仍然不对，应该是n-1后的数如果不是素数或其单素数次幂形式的数皆无解。最起码现在的帖子中没有反例。

mathe

If n ≡ 1 or 2 mod (4) and if the square free part of n contains at least one prime factor of the form 4k + 3 then there does not exist any projective plane of order n .

mathe

构造法

1. get_num_prim(p)=
   
2. {
   
3.     local(r,pass);
   
4.     r=factor(p-1);
   
5.     for(u=1,p-1,
   
6.         pass=u;
   
7.         for(v=1,length(r[,1]),
   
8.              if(Mod(u,p)^((p-1)/r[v,1])==1, pass=0;break())
   
9.         );
   
10.        if(pass > 0, break())
    
11.     );
    
12.     Mod(pass,p)
    
13. }
    
14. 
    
15. get_prim(p,n)=
    
16. {
    
17.     if(n==1,
    
18.         get_num_prim(p),
    
19.         ffgen(ffinit(p,n),a)
    
20.     )
    
21. }
    
22. 
    
23. irriducible(f)=
    
24. {
    
25.     local(r);
    
26.     r=factor(f);
    
27.     length(r[,1])==1 && r[1,2]==1
    
28. }
    
29. 
    
30. isprimitive(f,q)=
    
31. {
    
32.    local(r,pass);
    
33.    r=factor(q^3-1);pass=1;
    
34.    for(u=1,length(r[,1]),
    
35.         if(Mod(x,f)^((q^3-1)/r[u,1])==1,pass=0;break())
    
36.    );
    
37.    pass
    
38. }
    
39. 
    
40. findprimitive(a,p,n)=
    
41. {
    
42.     local(q,c1,c2,c3,f,found,b);
    
43.     q=p^n;c1=c2=0;c3=a/a;f=0;
    
44.     b=if(n==1,a,ffprimroot(a));
    
45.     for(u=1,q,
    
46.         c2=0;
    
47.         for(v=1,q,
    
48.             c3=a/a;
    
49.             for(w=1,q-1,
    
50.                 f=x^3+c1*x^2+c2*x+c3;
    
51.                 if(irriducible(f)&&isprimitive(f,q),
    
52.                     found=1;break()
    
53.                 );
    
54.                 c3=c3*b
    
55.             );
    
56.             if(found,break());
    
57.             if(v==1,c2=a/a,c2=c2*b);
    
58.         );
    
59.         if(found,break());
    
60.         if(u==1,c1=a/a,c1=c1*b);
    
61.     );
    
62.     f
    
63. }
    
64. 
    
65. findany(p,n)=
    
66. {
    
67.    local(f,q,v);
    
68.    a=get_prim(p,n);q=p^n;
    
69.    f=findprimitive(a,p,n);v=vector(0);
    
70.    if(f==0,print("Fail to find primitive function"),
    
71.        for(u=1, q^2+q,
    
72.              if(polcoeff(component(Mod(x,f)^u,2),2)==0,v=concat(v,u))
    
73.        )
    
74.    );
    
75.   kill(a);
    
76.   v
    
77. }
    
78. 
    
79. vecmatch(v1,v2)=
    
80. {
    
81.     local(matched);
    
82.     matched=1;
    
83.     for(u=1,length(v1),
    
84.          if(v1[u]<>v2[u], matched=0;break())
    
85.     );
    
86.     matched
    
87. }
    
88. 
    
89. vecrevert(w,m)=
    
90. {
    
91.      local(r);
    
92.      r=vector(length(w));
    
93.      r[1]=1;
    
94.      for(u=2,length(w),
    
95.           r[u]=m+1-w[length(w)+2-u]
    
96.      );
    
97.     r
    
98. }
    
99. 
    
100. vecshift(w,m)=
     
101. {
     
102.     local(h,r);h=length(w);
     
103.     if(w[1]==1,
     
104.           w,
     
105.           for(u=1,length(w)-1,
     
106.                if(w[u]+1==w[u+1],
     
107.                     h=u;break()
     
108.                )
     
109.           );
     
110.           r=vector(length(w));
     
111.           for(u=h+1,length(w),
     
112.               r[u-h]=(w[u]-w[h])%m
     
113.           );
     
114.           r[length(w)+1-h]=(-w[h])%m;
     
115.           for(u=1,h-1,
     
116.               r[length(w)+1-h+u]=(w[u]-w[h])%m
     
117.           );
     
118.          r
     
119.      )
     
120. }
     
121. 
     
122. findall(p,n)=
     
123. {
     
124.     local(v,q,m,r,w,rw,matched);
     
125.     q=p^n;m=q^2+q+1;
     
126.     r=matrix(0,q);
     
127.     v=findany(p,n);
     
128.     r=concat(r,v);
     
129.     w=vector(q);
     
130.     for(u=1,m-1,
     
131.          if(gcd(m,u)==1,
     
132.              for(k=1,q,
     
133.                  w[k]=(v[k]*u)%m;
     
134.              );
     
135.              w=vecsort(w);
     
136.              w=vecshift(w,m);
     
137.              rw=vecrevert(w,m);
     
138.              matched=0;
     
139.              for(k=1,length(r[,1]),
     
140.                   if(vecmatch(r[k,],w)||
     
141.                        vecmatch(r[k,],rw),
     
142.                       matched=1;break();
     
143.                   )
     
144.              );
     
145.              if(!matched, r=concat(r,w))
     
146.          )
     
147.     );
     
148.     r
     
149. }
     
150. 
     
151. normvec(v,m)=
     
152. {
     
153.     local(r);
     
154.     r=vector(length(v)+1);
     
155.     r[1]=1;
     
156.     for(u=1,length(v)-1,
     
157.          r[u+1]=v[u+1]-v[u]
     
158.     );
     
159.    r[length(v)+1]=m-v[length(v)];
     
160.    r
     
161. }
     
162. 
     
163. normresult(N)=
     
164. {
     
165.     local(p,n,q,r,m,rr);
     
166.     q=N-1;m=q^2+q+1;
     
167.     r=factor(q);
     
168.     if(length(r[,1])>1,
     
169.        [],
     
170.        p=r[1,1];n=r[1,2];
     
171.        r=findall(p,n);
     
172.        rr=matrix(length(r[,1]),q+1);
     
173.        for(u=1,length(r[,1]),
     
174.             rr[u,]=normvec(r[u,],m)
     
175.        );
     
176.        rr
     
177.    )
     
178. }

复制代码

如
(11:34) gp > normresult(8)
%216 =
[1 2 10 19  4  7  9  5]

[1 8 11  3 18 10  2  4]

[1 3  8  2 16  7 15  5]

[1 6 17 12  2 11  5  3]

[1 4 22  7  3  6  2 12]

[1 6 12  4 21  3  2  8]

mathe

wayne 发表于 2019-2-13 16:21
mathe能打破纪录吗，

经过数月等待，添加的3个零终于被OEIS接受了

chyanog

mathe已经完美解决了这个问题，还是发一下代码吧，不和C代码比，比前面的Mathematica代码要快得多，也是穷举，但是写法比较适合编译器优化，在6核笔记本上，计算n=8、9、10时都在1秒以内，n=12需要两分钟多。用了动态生成代码，预处理和编译也要花一点时间，建议搭配GCC，VC编译速度比较慢。

1. n=10;
   
2. max=n(n-1)/2+1
   
3. sum=n(n-1)+1
   
4. F=Symbol@FromCharacterCode[#+63487]&;
   
5. syms=Array[F,n-2];
   
6. syms=Join[{1},syms,{sum-1-Total@syms}]
   
7. vars=Union@@Table[Total/@Partition[syms,i,1,1],{i,n}];
   
8. cond=2<=Last@syms<=max&&And@@Thread[1<=Complement[vars,syms]<=sum]&&Simplify@LogicalExpand[Unequal@@vars];//AbsoluteTiming
   
9. 
   
10. iter=Table[{F[x+1],2,If[Select[cond,F[x]==Last@Sort@Cases[#,_Symbol,-1]&],Evaluate@If[x<n-2,max,2],0]},{x,n-2}];
    
11. 
    
12. cf=Compile[{{\[FormalA],_Integer}},
    
13. Module[{B=Internal`Bag[Rest@{0}]},
    
14. Do[Internal`StuffBag[B,#,1],##2];
    
15. Internal`BagPart[B,All]
    
16. ],CompilationTarget->"C",RuntimeOptions->"Speed",RuntimeAttributes->{Listable}
    
17. ]&[syms,Sequence@@iter];//AbsoluteTiming
    
18. 
    
19. Join@@(Partition[#,n]&/@cf[Range[2,max]])//AbsoluteTiming

复制代码

n=10时的结果