OFFSET
1,2
COMMENTS
It appears that starting at 64, the sequence follows the pattern, 2^(2n) + {0, 16, 21, 36, 64, 66, 84, 96, 144, 168} * 2^(2n-6) for n >= 3. - T. D. Noe, Jun 20 2012
LINKS
Robert Israel, Table of n, a(n) for n = 1..113 (first 100 terms from T. D. Noe)
FORMULA
Empirical g.f.: x*(1 +2*x +4*x^2 +5*x^3 +8*x^4 +10*x^5 +13*x^6 +16*x^7 +20*x^8 +25*x^9 +28*x^10 +29*x^11 +24*x^12 +32*x^13 +26*x^14 +24*x^15 +28*x^16 +21*x^17 +20*x^18 +28*x^19 +2*x^20) / ((1 -2*x^5)*(1 +2*x^5)). - Colin Barker, Feb 09 2016
MAPLE
# this requires Maple 17 or later
m:= 5000: # to find all entries <= m^2
N:= m^2:
with(SignalProcessing):
with(ArrayTools):
A1:= Array(0..N, datatype=float[8]):
for i from 1 to m do A1[i^2]:= 1 end do:
B:= Convolution(A1, A1):
A2:= Array(0..N, datatype=float[8]):
Copy(N+1, B, A2):
All1:= Array(0..N, datatype=float[8], fill=1):
MinimumEvery(A2, All1, container=A2):
B:= Convolution(A1, A2):
A3:= Array(0..N, datatype=float[8]):
Copy(N+1, B, A3);
MinimumEvery(A3, All1, container=A3);
MaximumEvery(A1, A2, container=A2);
B:= evalhf(map(round, A3-A2));
MinimumEvery(B, Array(0..N, datatype=float[8]), container=B);
R:= ArrayTools[SearchArray](B);
A000549:= map(`-`, convert(R, list), 1); # Robert Israel, Mar 28 2014
MATHEMATICA
A000549 = Reap[For[n=1, n <= 10^5, n++, If[SquaresR[2, n] != 0, If[Union[ PowersRepresentations[n, 3, 2][[All, 1]]] == {0}, Print[n]; Sow[n]]]] ][[2, 1]] (* Jean-François Alcover, Feb 09 2016 *)
PROG
(MATLAB) Asq = zeros(1, 10^7);
Asq( [0:3162] .^ 2 + 1) = 1;
Psq = zeros(1, 10^7);
Psq( [1:3162] .^ 2 + 1) = 1;
As2 = conv(Asq, Asq);
As2 = min(1, As2(1:10^7));
Ps2 = conv(Psq, Psq);
Ps2 = min(1, Ps2(1:10^7));
As3 = conv(Asq, As2);
As3 = min(1, As3(1:10^7));
Ps3 = conv(Psq, Ps2);
Ps3 = min(1, Ps3(1:10^7));
D = As3 - Ps3;
A000549 = find(D(2:end))
% Robert Israel, Mar 26 2014
(PARI) istwo(n)=if(n<3, return(n>=0)); my(f=factor(n>>valuation(n, 2))); for(i=1, #f~, if(f[i, 2]%2&&f[i, 1]%4==3, return(0))); 1
is(n)=if(!istwo(n), return(0)); if(n<3, return(1)); for(x=sqrtint(n\3), sqrtint(n-2), my(t=n-x^2); for(y=sqrtint(t\2), sqrtint(t-1), if(t&&issquare(t-y^2), return(0)))); 1 \\ Charles R Greathouse IV, Jan 28 2016
CROSSREFS
KEYWORD
nonn
AUTHOR
EXTENSIONS
Extended by T. D. Noe, Jun 20 2012
STATUS
approved