View Single Post
Old 2022-04-17, 18:03   #2
mathwiz
 
Mar 2019

1001101102 Posts
Default

How about:

Code:
?  forprime(p = 2, 10000000, q=p*p+p+1;  if(isprime(q) && vecmax(factorint(q^3-1)[, 1]) == p, print(p)))
125669
138209
254537
309629
532187
1107497
1126523
1210103
1225817
1287329
1524431
1534349
1539719
1720181
1793123
1814609
1861151
1920731
1932071
1974881
2270423
2366057
2490479
2494931
2530373
2586377
2725841
2755943
2782667
2885837
...
mathwiz is offline   Reply With Quote