OFFSET
1,2
COMMENTS
A number k has a coprime divisor shift s if GCD(d + s, k) = 1 for all divisors d of k.
If k is in the sequence, then all divisors of k are in the sequence too.
From David A. Corneth, Oct 06 2023: (Start)
As a consequence of the above, if some m is not in the sequence then any multiple of m is not in the sequence either. So all terms are odd as 2 is not in the sequence.
To see if k is a term we can do the following:
- Check whether k is even; if so then k is not in the sequence.
- Make a list of the divisors of k. Let tau(k) be the number of divisors of k. Then for each prime p <= tau(k) if the number of residues mod p of the divisors of k is equal to p then k is not in the sequence. Otherwise by the Chinese Remainder Theorem we can find a number s such that gcd(d + s, n) = 1 for all d.
So every odd prime is a term. (End)
LINKS
M. Farrokhi D. G., Table of n, a(n) for n = 1..10000
EXAMPLE
1 is a term since GCD(1 + 0, 1) = 1.
2 is not a term since GCD(1 + s, 2) > 1 or GCD(2 + s, 2) > 1 for all nonnegative integers s.
3 is a term since GCD(1 + 1, 3) = GCD(3 + 1, 3) = 1.
MAPLE
aList := proc(len) local isds, findds;
isds := (k, s) -> andmap(d -> igcd(d + s, k) = 1, NumberTheory:-Divisors(k));
findds := k -> ormap(s -> isds(k, s), [seq(1..k)]);
select(k -> findds(k), [seq(1..len)]) end:
aList(133); # Peter Luschny, Oct 06 2023
# More efficient, after David A. Corneth:
isa := proc(n) local d, p, t; p := 3;
if irem(n, 2) = 0 then return false fi;
d := NumberTheory:-Divisors(n);
while p < nops(d) do
{seq(irem(t, p), t = d)};
if nops(%) = p then return false fi;
p := nextprime(p);
od: true end:
aList := len -> select(isa, [seq(1..len)]):
aList(145); # Peter Luschny, Oct 07 2023
MATHEMATICA
isa[n_] := Module[{d, p, t}, p = 3; If[Mod[n, 2] == 0, Return[False]]; d = Divisors[n]; While[p < Length[d], If[Length[Union@Table[Mod[t, p], {t, d}]] == p, Return[False]]; p = NextPrime[p]]; True];
aList[len_] := Select[Range[len], isa];
aList[145] (* Jean-François Alcover, Oct 27 2023, after Peter Luschny *)
PROG
(PARI)
isds(k, s)={fordiv(k, d, if(gcd(d+s, k)<>1, return(0))); 1}
findds(k)={for(s=1, k, if(isds(k, s), return(s))); 0}
select(k->findds(k), [1..150]) \\ Andrew Howroyd, Oct 05 2023
(PARI)
is(n) = {
if(!bitand(n, 1), return(0));
my(d = divisors(n));
forprime(p = 3, #d,
if(#Set(d % p) == p,
return(0)
)
); 1
} \\ David A. Corneth, Oct 06 2023
(GAP)
CoprimeDivisorShift := function(n)
local shift, divisors;
shift := 0;
if not IsPrimeInt(n) and First(PrimeDivisors(n), p -> CoprimeDivisorShift(n / p) = infinity) <> fail then
shift := infinity;
fi;
if shift < infinity then
divisors := DivisorsInt(n);
while shift < n and First(divisors, d -> GcdInt(d + shift, n) > 1) <> fail do
shift := shift + 1;
od;
if shift = n then
shift := infinity;
fi;
fi;
return shift;
end;
CROSSREFS
KEYWORD
nonn,easy
AUTHOR
M. Farrokhi D. G., Oct 05 2023
STATUS
approved