[seqfan] reducing Devaraj non-Carmichel numbers in A164946

[Richard Mathar]

I have problems reproducing the numbers in http://oeis.org/A164946. I get
73947,
250811,
604803,
2237805,
2544723,
35839470,
24413481,
8039307,
8213956,
95828168,
107785924,
39242222949303,
12827318,
114397938,
14205361,
19722131,
81740112592896,
193098816,
26308112,
215405768,
125012678906250,
136759472331250,
236301822947449,
269517889,
287152344,
157098832,

if we start from the list in A104017, then take the recipe in A104016
(there is no more precise description avail as it seems, because
A164946 says "computation... is similar...."): compute the product
of the square of the greatest common divisor of the 1-reduced prime factors,
multiply by the power of (n-1) and divide through the product of the 1-reduced
prime factors. The current entries are usually (but not always) a factor 4
larger than the results above.

First example is 11305 = 5*7*17*19, prime-factors are 4,6,16,18 after
reduction, with gcd(4,6,16,18) = 2. So 2^2*(11305-1)^2/4/6/16/18 = 73947
according to my interpretation.

In Maple:

DevRedu := proc(n)
	local d,dS,a,l;

	# list of prime factors
	d := convert(numtheory[factorset](n),list) ;

	# the variable r in A104016 ;
	dS := nops(d) ;

	# square of gcd(...)^2*(N-1)^(r-2) as in A104016
	a := igcd(seq(l-1,l=d))^2*(n-1)^(dS-2) ;

	# divide through (p1-1)*...(pr-1) as in A104016
	for l in d do
		a := a/(l-1) ;
	end do:
end proc:
L :=[
11305, 39865, 96985, 401401, 464185, 786961, 1106785, 1296505, 1719601, 1993537, 2242513, 2615977, 2649361, 2722681, 3165961, 3181465, 3755521, 4168801, 4229601, 4483297, 4698001, 5034601, 5381265, 5910121, 5977153, 7177105
] ;
for l in L do
	printf("%d,\n",DevRedu(l)) ;
end do: