
PROG

(PARI)
/* define $cy(m, n) = Phi_m(n)$ the $m$th cyclotomic polynomial evaluated at $t=n$ */
cy(m, n) = {local(po); po = polcyclo(m, t); subst(po, t, n); }
/* for fixed $m$ compute $cy(m, n)$ */
cy1(n) = {subst(po, t, n); }
/* search from m <a to b possible solutions m, n, using parity of m */
scy1(a, b, fr) =
{local(m, n, r, c, c1, c2, d, g, po, be, t1, t2, t3, st, dd);
st = 1; dd = 1; be = ba; t1 = gettime();
for(m=a, b, po = polcyclo(m, t); if(m % 2 == 0, st = 2; dd = 3; ,
st = 1; dd = 2; );
forstep(n=dd, m1, st, g = gcd(n, m);
if(g == 1, c1 = cy1(n); r = 2*m*n+1;
if(c1 % r == 0, c2 = cy(n, m); d = gcd(c2, c1);
if(d == r, print([m, n, r]);
); ); ); );
if(m % fr == 0,
t2 = gettime(); t3 = t2+t1; print([m, t3, ((ma)/be)*100.0]);
t1 = t3;
); ); }
