OFFSET
1,3
COMMENTS
a(n) = 1 if n is a term in A005574 (numbers n such that n^2 + 1 is prime).
a(n) = 1 if gpf(k^2 + 1) <> gpf(n^2 + 1) for every positive integer k < n.
LINKS
Michel Lagneau, Table of n, a(n) for n = 1..10000
EXAMPLE
a(3) = 2 because the greatest prime divisor of 3^2 + 1 is 5 and n=3 is the 2nd positive value of n at which gpf(n^2 + 1) = 5; the 1st is n=2: gpf(2^2 + 1) = 5.
a(313) = 7 because the greatest prime divisor of 313^2 + 1 is 101, and n=313 is the 7th positive value of n at which this occurs:
10^2 + 1 = 101;
91^2 + 1 = 2 * 41 * 101;
111^2 + 1 = 2 * 61 * 101;
192^2 + 1 = 5 * 73 * 101;
212^2 + 1 = 5 * 89 * 101;
293^2 + 1 = 2 * 5^2 * 17 * 101;
313^2 + 1 = 2 * 5 * 97 * 101.
MAPLE
with(numtheory):nn:=200:T:=array(1..nn):k:=0:
for m from 1 to nn do:
x:=factorset(m^2+1):n1:=nops(x):p:=x[n1]:k:=k+1:T[k]:=p:
od:
for n from 1 to 150 do:
q:=T[n]:ii:=0:
for i from 1 to n do:
if T[i]=q then ii:=ii+1:
else
fi:
od:
printf(`%d, `, ii):
od:
# Simpler version:
N:= 1000: # to get a(n) for n <= N
T:= Array(1..N):
for n from 1 to N do
T[n]:= max(numtheory:-factorset(n^2+1));
A[n]:= numboccur(T, T[n]);
od:
seq(A[n], n=1..N); # Robert Israel, Aug 12 2014
PROG
(PARI) a(n) = my(gn = vecmax(factor(n^2+1)[, 1])); sum(k=1, n, vecmax(factor(k^2+1)[, 1]) == gn); \\ Michel Marcus, Sep 10 2017
CROSSREFS
KEYWORD
nonn
AUTHOR
Michel Lagneau, Aug 11 2014
EXTENSIONS
Edited by Jon E. Schoenfield, Sep 10 2017
STATUS
approved