? for(q=1,1000,for(n=1,1000,m=q-n^2;if(m>0,s=(7^q+1)/(7^m+1);if(denominator(s)==
1,if(ispseudoprime(s)==1,print([m,q,n]))))))
[1, 17, 4]
[8, 24, 4]
[2, 38, 6]
[4, 148, 12]

? z=prod(j=1,3*10^4,prime(j));for(q=1,10^4,for(n=1,10^4,m=q-n^2;if(m>0,s=(7^q+1)
/(7^m+1);if(denominator(s)==1,if(gcd(s,z)==1,if(ispseudoprime(s)==1,print([m,q,n
])))))))
[1, 17, 4]
[8, 24, 4]
[2, 38, 6]
[4, 148, 12]
?