OFFSET
1,2
LINKS
FORMULA
MAPLE
isA000404 := proc(n)
local x, y ;
for x from 1 do
if x^2> n then
return false;
end if;
for y from 1 do
if x^2+y^2 > n then
break;
elif x^2+y^2 = n then
return true;
end if;
end do:
end do:
end proc:
A018825 := proc(n)
if n = 1 then
1;
else
for a from procname(n-1)+1 do
if not isA000404(a) then
return a;
end if;
end do:
end if;
end proc:
seq(A018825(n), n=1..30) ; # R. J. Mathar, Jul 28 2014
MATHEMATICA
q=13; q2=q^2+1; lst={}; Do[Do[z=a^2+b^2; If[z<=q2, AppendTo[lst, z]], {b, a, 1, -1}], {a, q}]; lst; u=Union@lst; Complement[Range[q^2], u] (* Vladimir Joseph Stephan Orlovsky, May 30 2010 *)
PROG
(Haskell)
import Data.List (elemIndices)
a018825 n = a018825_list !! (n-1)
a018825_list = tail $ elemIndices 0 a025426_list
-- Reinhard Zumkeller, Aug 16 2011
(PARI) is(n)=my(f=factor(n), t=prod(i=1, #f~, if(f[i, 1]%4==1, f[i, 2]+1, if(f[i, 2]%2 && f[i, 1]>2, 0, 1)))); if(t!=1, return(!t)); for(k=sqrtint((n-1)\2)+1, sqrtint(n-1), if(issquare(n-k^2), return(0))); 1 \\ Charles R Greathouse IV, Sep 02 2015
CROSSREFS
KEYWORD
nonn
AUTHOR
STATUS
approved