|
%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_
|
