|
PROG
|
(PARI)
{
\\Program runs for arbitrary B.\\
B=10^12; N=primepi(B/(627*17*19*23));
p=vector(N, n, prime(n));
in=primepi((B/210)^(1/3));
P=prod(i=1, 27, p[i]); Q=prod(i=28, in, p[i]);
v=28; d=[[1, 2], [-1, 1], [-2, -1]]; i3=6;
while(6*p[i3]^5<B,
k1=p[i3];
for(i2=2, i3-1,
k1*=p[i2];
for(i1=1, i2-1,
k1*=p[i1];
if(k1>626,
if(k1*p[i3+1]*p[i3+2]*p[i3+3]*p[i3+4]<B,
if(6*p[v-1]<k1,
while(6*p[v]<k1, Q/=p[v]; P*=p[v]; v++),
while(6*p[v-1]>=k1, v--; Q*=p[v]; P/=p[v]));
r=(B\k1)^(1/4); j1=i3+1;
while(p[j1]<r,
k2=p[j1]*k1; s(B\k2)^(1/3); j2=j1+1;
while(p[j2]<s,
k3=p[j2]*k2; t=sqrt(B\k3); j3=j2+1;
while(p[j3]<t,
k4=k3*p[j3]; j4=j3+1; u=B\k4;
while(p[j4]<=u,
K=k4*p[j4];
a=K+d[K%4][1]; b=K+d[K%4][2];
a1=gcd(P, a); y1=omega(a1);
if(y1>2,
f=1; if(y1==3, if(a1>j1, f=0));
if(f,
b1=gcd(P, b); z1=omega(b1);
if(z1>2,
if(z1==3, if(b1>j1, f=0));
if(f,
a2=a/a1;
if(gcd(a1, a2)==1,
b2=b/b1;
if(gcd(b1, b2)==1,
a21=gcd(a2, Q); a22=a2/a21;
if(gcd(a21, a22)==1,
y=y1+omega(a21);
if(y>4,
if(y<8,
b21=gcd(Q, b2); b22=b2/b21;
if(gcd(b21, b22)==1,
z=z1+omega(b21);
if(z>4,
if(z<8,
if(y+omega(a22)==7,
if(z+omega(b22)==7,
f1=factor(a1);
if(f1[1, 1]*f1[2, 1]*f1[3, 1]<k1,
f2=factor(b1);
if(f2[1, 1]*f2[2, 1]*f2[3, 1]<k1,
print1((K\4)*4+2" ")
))))))))))))))));
j4++);
j3++);
j2++);
j1++)));
k1/=p[i1]);
k1/=p[i2]);
i3++)
}
|