A268629 Primes p that have no squareful primitive roots less than p.

3, 5, 7, 13, 17, 19, 23, 31, 41, 43, 47, 61, 71, 73, 79, 97, 103, 127, 191, 193, 223, 239, 241, 311, 313, 337, 409, 433, 439, 457, 479, 601, 719, 769, 839, 911, 1009, 1031, 1033, 1129, 1151, 1201, 1249, 1319, 1321, 1559, 1801, 2089, 2281, 2521, 2689, 2999, 3049, 3361, 3529, 3889

Offset: 1

Author

Michel Marcus, Feb 09 2016

Keywords

• nonn

Examples

			The primitive roots of 7 less than 7 are 3 and 5. None of them are squareful so 7 is in the sequence.
8 is a primitive root of 11, and 8 is squareful, so 11 is not in the sequence.
		

Links

• Robert Israel, Table of n, a(n) for n = 1..114
• Stephen D. Cohen and Tim Trudgian, On the least square-free primitive root modulo p, arXiv:1602.02440 [math.NT], 2016.

Crossrefs

Cf. A001694, A001918, A060749.

Programs

Maple

N:= 10^6: # for terms <= N
S:= {1}: p:= 1:
do
  p:= nextprime(p);
  if p^2 > N then break fi;
  S:= S union map(t -> seq(t*p^i, i=2..floor(log[p](N/t))), select(`<=`,S,N/p^2));
od:
S:= sort(convert(S,list)):
nS:= nops(S):
filter:= proc(p) local i;
  if not isprime(p) then return false fi;
  for i from 1 to nS while S[i] < p do
    if numtheory:-order(S[i],p) = p-1 then return false fi
  od;
  true
end proc:
select(filter, [seq(i,i=3..N,2)]); # Robert Israel, Oct 27 2020

Mathematica

selQ[p_] := NoneTrue[PrimitiveRootList[p], #
= 2&]&];
Select[Prime[Range[2, 500]], selQ] (* Jean-François Alcover, Sep 28 2018 *)

PARI

ar(p) = my(r, pr, j); r=vector(eulerphi(p-1)); pr=znprimroot(p); for(i=1, p-1, if(gcd(i, p-1)==1, r[j++]=lift(pr^i))); vecsort(r) ; \\ from A060749
isok(p) = {my(v = ar(p)); for (i=1, #v, if (ispowerful(v[i]), return(0));); 1;}
lista(nn) = forprime(p=1, nn, if (isok(p), print1(p, ", ")));

Python

from functools import cache
from math import gcd
from itertools import count, islice
from sympy import factorint, prime, n_order
@cache
def is_squareful(n): return n == 1 or min(factorint(n).values()) > 1
def A268629_gen(): # generator of terms
    for n in count(1):
        p = prime(n)
        for i in range(1,p):
            if gcd(i,p) == 1 and is_squareful(i) and n_order(i, p)==p-1:
                break
        else:
            yield p
A268629_list = list(islice(A268629_gen(),20)) # Chai Wah Wu, Sep 14 2022