У меня тогда два варианта кода получилось для сбора-проверки подходящих чисел. Один на основе массива Set; поиск в нём быстр, но пополнение его очень не эффективно, т.к. каждый раз пересорт массива происходит.
Код:
pkh()=
{
local(ii, n, nl, p, pmax, k, h2, h3, g, l, z1, z2, o, dz, a, m, mm, a2, m2):int;
local(s):bool;
local(v, w, Y, Z):vec;
local(pr, st);
Y= read("pkh.dbt"):vec;
st= Set(Y[,1]);
Z= read("pkhset.dbt"):vec;
\\st= Set():vec;
\\Y= matrix(0,3):vec;
\\n= (5):int;
\\n= (52*10^5+1):int;
\\n= (10527*10^5+1):int;
\\n= (47907*10^5+1):int;
\\n= (52637*10^5+1):int;
\\n= (74121*10^5+1):int;
\\n= (407433*10^5+1):int;
\\n= (102116*10^6+1):int;
\\n= (125641*10^6+1):int;
n= (10^5+1):int;
\\pmax= 0:int;
pmax= vecmax(Y[,1]):int;
while(n < (10^6),
if(n%3, if(n%7, if(n%11, if(n%13, if(n%17, if(n%31, if(n%37,
if(issquarefree(n),
s= 1:bool;
w= factorint(n):vec;
nl= #w[,1]:int;
for(j=1, nl,
if(w[j,1]<=pmax,
if(!setsearch(st,w[j,1]), s= 0:bool; break)
)
);
if(s,
for(j=1, nl,
p= w[j,1]:int;
if(p<=pmax,
ii= setsearch(st,p):int;
a2= Z[ii,2]:int; m2= Z[ii,3]:int
,
g= 0:int;
h2= znorder(Mod(2,p)):int; h3= znorder(Mod(3,p)):int;
if(h2<>h3,
if(h2 % h3, s= 0:bool; break,
g= (h2/h3):int;
if(!issquarefree(g), s= 0:bool; break,
v= factorint(g):vec;
l= #v[,1]:int;
for(i=1, l,
if(!setsearch(st,v[i,1]), s= 0:bool; break(2))
)
)
)
);
pr= znprimroot(p);
z1= znlog(3:int, pr):int; z2= znlog(2:int, pr):int;
o= eulerphi(p):int; dz= gcd(o,z2):int;
if(z1 % dz, s= 0:bool; break(), k= ((z1/z2) % (o/dz)):int);
if(!(h2%24),
if(kronecker(6,lift(Mod(k,24)))==-1, s= 0:bool; print(n); break)
);
a2= (p*k):int; m2= (p*h2):int;
if(g>1,
for(i=1, l,
ii= setsearch(st,v[i,1]):int;
a= Z[ii,2]:int; m= Z[ii,3]:int;
mm= gcd(m,m2);
if(Mod(a,mm) == Mod(a2,mm),
a2= lift(chinese(Mod(a,m),Mod(a2,m2))):int;
m2= (m*m2/mm):int, s= 0:bool; break(2)
)
)
);
Y= concat(Y, [p, a2, m2]):vec;
pmax= p:int;
);
if(j==1, a= a2:int; m= m2:int,
mm= gcd(m,m2);
if(Mod(a,mm) == Mod(a2,mm),
a= lift(chinese(Mod(a,m),Mod(a2,m2))):int;
m= (m*m2/mm):int, s= 0:bool; break
)
);
)
);
if(s,
if(a>10^5,
if(a%2, if(a%3, if(a%7, if(a%11, if(a%13, if(a%17, if(a%31, if(a%37,
if(lift(Mod(2,a)^a)==3,
print("a= ",a," ",factorint(a)," n= ",n," ",factorint(n));
write("pkh.txt",a," ",factorint(a)," ",n," ",factorint(n))
)
))))))))
)
)
)
)))))));
n+=2
);
write("pkhnew.dbt", Y);
}
pkhset()=
{
local(l):int;
local(Y, Z, st):vec;
Y= read("pkh.dbt"):vec;
st= Set(Y[,1]):vec;
l= (#st):int;
Z= matrix(l,3):vec;
for(i=1, l, for(j=1, l, if(eval(st[i])==Y[j,1], Z[i,]= Y[j,]; print(i); break())));
write("pkhset.dbt", Z);
}
Другой на основе собственного бинарного поиска в массиве.
Код:
zsearch(p:int, Y:vec)=
{
local(j, i1, i2, it, t):int;
if(p<5, j= 0:int,
i1= 1:int;
i2= (#Y):int;
t= -1:int;
while(p!=t,
it= ((i2+i1)\2):int;
t= Y[it]:int;
if(p<t, i2= (it-1):int,
if(p>t, i1= (it+1):int,
j= it:int; break
)
);
if((i1>i2), j= 0:int; break)
)
);
return(j:int)
};
pkh()=
{
local(Z, w):vec;
local(n, g, i, h, h3, z2, z3, o, dz, k, a, m, mm, a1, m1, a2, m2, d, L):int;
local(s):bool;
local(pr);
\\Z= matrix(0, 3):vec;
Z= read("mz78.dbt"):vec;
L= (#Z~):int;
if(L==0, n= (5):int, n= (Z[L,1]+2):int);
\\n= (526*10^4+1):int;
\\n= (105274*10^4+1):int;
\\n= (479070*10^4+1):int;
\\n= (741213*10^4+1):int;
while(n < 10^8,
if(n%3, if(n%7, if(n%11, if(n%13, if(n%17, if(n%31, if(n%37,
if(issquarefree(n),
s= 1:bool;
if(isprime(n),
g= 0:int;
h= znorder(Mod(2, n)):int; h3= znorder(Mod(3, n)):int;
if(h<>h3,
if(h%h3, s= 0:bool,
g= (h/h3):int;
i= zsearch(g, Z[,1]):int;
if(!i, s= 0:bool)
)
);
if(s,
pr= znprimroot(n);
z2= znlog(2:int, pr):int; z3= znlog(3:int, pr):int;
o= eulerphi(n):int; dz= gcd(o, z2):int;
if(z3%dz, s= 0:bool,
k= ((z3/z2)%(o/dz)):int;
if(!(h%24), if(kronecker(6, lift(Mod(k, 24)))==-1, s= 0:bool; print(n)))
)
);
if(s,
a= (n*k):int; m= (n*h):int;
if(g,
a2= (g*Z[i,2]):int; m2= (g*Z[i,3]):int;
mm= gcd(m, m2):int;
if(Mod(a, mm) == Mod(a2, mm),
a= lift(chinese(Mod(a, m), Mod(a2, m2))):int;
m= (m*m2/mm):int, s= 0:bool
)
)
)
,
w= divisors(n):vec;
d= w[2]:int;
i= zsearch(d, Z[,1]):int;
if(i, a1= (d*Z[i,2]):int; m1= (d*Z[i,3]):int, s= 0:bool);
if(s,
L= (#w):int;
d= w[L-1]:int;
i= zsearch(d, Z[,1]):int;
if(i, a2= (d*Z[i,2]):int; m2= (d*Z[i,3]):int, s= 0:bool)
);
if(s,
mm= gcd(m1, m2):int;
if(Mod(a1, mm) == Mod(a2, mm),
a= lift(chinese(Mod(a1, m1), Mod(a2, m2))):int;
m= (m1*m2/mm):int, s= 0:bool
)
);
);
if(s,
if(a%2, if(a%3, if(a%7, if(a%11, if(a%13, if(a%17, if(a%31, if(a%37,
if(lift(Mod(2,a)^a)==3,
print("a= ",a," ",factorint(a)," n= ",n," ",factorint(n));
write("pkh.txt",a," ",factorint(a)," ",n," ",factorint(n))
)
))))))));
Z= concat(Z, [n,a/n,m/n])
);
)
)))))));
n+=2
);
write("mz8.dbt", Z)
};
Мне хотелось создать код, который бы всю базу подходящих чисел держал в памяти, быстро её пополнял и быстро в ней искал. Нужно было выбрать, с каким массивом лучше работать - Set, вектор или матрицу; соответственно в каком формате сохранять. Нужно минимизировать дисковые операции, лучше однократно при входе читать из файла базы и при выходе писать данные в файл базы. Возможно ли сделать управляемый выход, т.е. чтоб в любой момент вручную остановить исполнение и корректно сохранить набранные данные? Какие слабые места наблюдаются в этом коде? Если кто понял алгоритм
Руста, то не могли бы объяснить его на конкретном примере?