A018825 - OEIS

%I #29 Sep 02 2015 11:40:54

%S 1,3,4,6,7,9,11,12,14,15,16,19,21,22,23,24,27,28,30,31,33,35,36,38,39,

%T 42,43,44,46,47,48,49,51,54,55,56,57,59,60,62,63,64,66,67,69,70,71,75,

%U 76,77,78,79,81,83,84,86,87,88,91,92,93,94,95,96,99,102,103,105,107,108,110

%N Numbers that are not the sum of 2 nonzero squares.

%H T. D. Noe, <a href="/A018825/b018825.txt">Table of n, a(n) for n = 1..10000</a>

%H <a href="/index/Su#ssq">Index entries for sequences related to sums of squares</a>

%F A025426(a(n)) = 0; A063725(a(n)) = 0. - _Reinhard Zumkeller_, Aug 16 2011

%p isA000404 := proc(n)

%p     local x,y ;

%p     for x from 1 do

%p         if x^2> n then

%p             return false;

%p         end if;

%p         for y from 1 do

%p             if x^2+y^2 > n then

%p                 break;

%p             elif x^2+y^2 = n then

%p                 return true;

%p             end if;

%p         end do:

%p     end do:

%p end proc:

%p A018825 := proc(n)

%p     if n = 1 then

%p         1;

%p     else

%p         for a from procname(n-1)+1 do

%p             if not isA000404(a) then

%p                 return a;

%p             end if;

%p         end do:

%p     end if;

%p end proc:

%p seq(A018825(n),n=1..30) ; # _R. J. Mathar_, Jul 28 2014

%t 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 *)

%o (Haskell)

%o import Data.List (elemIndices)

%o a018825 n = a018825_list !! (n-1)

%o a018825_list = tail $ elemIndices 0 a025426_list

%o -- _Reinhard Zumkeller_, Aug 16 2011

%o (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

%Y Cf. A022544, A081324, A000404 (complement), A004431.

%K nonn

%O 1,2

%A _David W. Wilson_