This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.
composite(n) = my(i=0, c=2); while(1, if(!ispseudoprime(c), i++); if(i==n, return(c)); c++) a273339(n) = forcomposite(c=1, , if(Mod(n, c^2)^(c-1)!=1, return(c))) a(n) = my(b=2, c=composite(n)); while(a273339(b)!=c, b++); b
15 satisfies the congruence 26^(15-1) == 1 (mod 15^2) and 15 < 26, so 26 is a term of the sequence.
N:= 1000: # to get all terms <= N
Res:= {}:
for c from 4 to N-1 do
if not isprime(c) then
for m in map(rhs@op, [msolve(x^(c-1)-1, c^2)]) do
if m > c and m <= N then Res:= Res union {m, seq(k*c^2+m, k=1..(N-m)/c^2)}
else Res:= Res union {seq(k*c^2+m, k=1..(N-m)/c^2)}
fi
od
fi
od:
sort(convert(Res,list)); # Robert Israel, Apr 20 2017
nn = 361; c = Select[Range@ nn, CompositeQ]; Select[Range@ nn, Function[n, Count[TakeWhile[c, # <= n &], k_ /; Mod[n^(k - 1), k^2] == 1] > 0]] (* Michael De Vlieger, May 30 2016 *)
is(n) = forcomposite(c=1, n-1, if(Mod(n, c^2)^(c-1)==1, return(1))); return(0)
Comments