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 =
? factor(2^15-1) //因子分解$F_{r^3}$的阶
%85 =
[ 31 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后的数如果不是素数或其单素数次幂形式的数皆无解。最起码现在的帖子中没有反例。
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 .
构造法
get_num_prim(p)=
{
local(r,pass);
r=factor(p-1);
for(u=1,p-1,
pass=u;
for(v=1,length(r[,1]),
if(Mod(u,p)^((p-1)/r)==1, pass=0;break())
);
if(pass > 0, break())
);
Mod(pass,p)
}
get_prim(p,n)=
{
if(n==1,
get_num_prim(p),
ffgen(ffinit(p,n),a)
)
}
irriducible(f)=
{
local(r);
r=factor(f);
length(r[,1])==1 && r==1
}
isprimitive(f,q)=
{
local(r,pass);
r=factor(q^3-1);pass=1;
for(u=1,length(r[,1]),
if(Mod(x,f)^((q^3-1)/r)==1,pass=0;break())
);
pass
}
findprimitive(a,p,n)=
{
local(q,c1,c2,c3,f,found,b);
q=p^n;c1=c2=0;c3=a/a;f=0;
b=if(n==1,a,ffprimroot(a));
for(u=1,q,
c2=0;
for(v=1,q,
c3=a/a;
for(w=1,q-1,
f=x^3+c1*x^2+c2*x+c3;
if(irriducible(f)&&isprimitive(f,q),
found=1;break()
);
c3=c3*b
);
if(found,break());
if(v==1,c2=a/a,c2=c2*b);
);
if(found,break());
if(u==1,c1=a/a,c1=c1*b);
);
f
}
findany(p,n)=
{
local(f,q,v);
a=get_prim(p,n);q=p^n;
f=findprimitive(a,p,n);v=vector(0);
if(f==0,print("Fail to find primitive function"),
for(u=1, q^2+q,
if(polcoeff(component(Mod(x,f)^u,2),2)==0,v=concat(v,u))
)
);
kill(a);
v
}
vecmatch(v1,v2)=
{
local(matched);
matched=1;
for(u=1,length(v1),
if(v1<>v2, matched=0;break())
);
matched
}
vecrevert(w,m)=
{
local(r);
r=vector(length(w));
r=1;
for(u=2,length(w),
r=m+1-w
);
r
}
vecshift(w,m)=
{
local(h,r);h=length(w);
if(w==1,
w,
for(u=1,length(w)-1,
if(w+1==w,
h=u;break()
)
);
r=vector(length(w));
for(u=h+1,length(w),
r=(w-w)%m
);
r=(-w)%m;
for(u=1,h-1,
r=(w-w)%m
);
r
)
}
findall(p,n)=
{
local(v,q,m,r,w,rw,matched);
q=p^n;m=q^2+q+1;
r=matrix(0,q);
v=findany(p,n);
r=concat(r,v);
w=vector(q);
for(u=1,m-1,
if(gcd(m,u)==1,
for(k=1,q,
w=(v*u)%m;
);
w=vecsort(w);
w=vecshift(w,m);
rw=vecrevert(w,m);
matched=0;
for(k=1,length(r[,1]),
if(vecmatch(r,w)||
vecmatch(r,rw),
matched=1;break();
)
);
if(!matched, r=concat(r,w))
)
);
r
}
normvec(v,m)=
{
local(r);
r=vector(length(v)+1);
r=1;
for(u=1,length(v)-1,
r=v-v
);
r=m-v;
r
}
normresult(N)=
{
local(p,n,q,r,m,rr);
q=N-1;m=q^2+q+1;
r=factor(q);
if(length(r[,1])>1,
[],
p=r;n=r;
r=findall(p,n);
rr=matrix(length(r[,1]),q+1);
for(u=1,length(r[,1]),
rr=normvec(r,m)
);
rr
)
}
如
(11:34) gp > normresult(8)
%216 =
wayne 发表于 2019-2-13 16:21
mathe能打破纪录吗，
经过数月等待，添加的3个零终于被OEIS接受了
本帖最后由 chyanog 于 2019-12-25 21:35 编辑
mathe已经完美解决了这个问题，还是发一下代码吧，不和C代码比，比前面的Mathematica代码要快得多，也是穷举，但是写法比较适合编译器优化，在6核笔记本上，计算n=8、9、10时都在1秒以内，n=12需要两分钟多。用了动态生成代码，预处理和编译也要花一点时间，建议搭配GCC，VC编译速度比较慢。n=10;
max=n(n-1)/2+1
sum=n(n-1)+1
F=Symbol@FromCharacterCode[#+63487]&;
syms=Array;
syms=Join[{1},syms,{sum-1-Total@syms}]
vars=Union@@Table,{i,n}];
cond=2<=Last@syms<=max&&And@@Thread<=sum]&&Simplify@LogicalExpand;//AbsoluteTiming
iter=Table[{F,2,If==Last@Sort@Cases[#,_Symbol,-1]&],Evaluate@If,0]},{x,n-2}];
cf=Compile[{{\,_Integer}},
Module[{B=Internal`Bag},
Do,##2];
Internal`BagPart
],CompilationTarget->"C",RuntimeOptions->"Speed",RuntimeAttributes->{Listable}
]&;//AbsoluteTiming
Join@@(Partition[#,n]&/@cf])//AbsoluteTiming
n=10时的结果